Paper figures
In [1]:
Copied!
import pandas as pd
import matplotlib
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
import seaborn as sns
from datetime import datetime
import numpy as np
from pathlib import Path
from matplotlib.colors import LightSource
import pygmt
# sns.set_theme()
# %matplotlib widget
params = {
'font.size': 18,
'axes.labelsize': 18,
'axes.titlesize': 18,
'legend.fontsize': 18,
'xtick.labelsize': 18,
'ytick.labelsize': 18,
'image.origin': "lower",
"figure.figsize": (4 * 2.5, 3 * 2.5),
"savefig.bbox": "tight",
"savefig.dpi": 300,}
matplotlib.rcParams.update(params)
import pandas as pd
import matplotlib
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
import seaborn as sns
from datetime import datetime
import numpy as np
from pathlib import Path
from matplotlib.colors import LightSource
import pygmt
# sns.set_theme()
# %matplotlib widget
params = {
'font.size': 18,
'axes.labelsize': 18,
'axes.titlesize': 18,
'legend.fontsize': 18,
'xtick.labelsize': 18,
'ytick.labelsize': 18,
'image.origin': "lower",
"figure.figsize": (4 * 2.5, 3 * 2.5),
"savefig.bbox": "tight",
"savefig.dpi": 300,}
matplotlib.rcParams.update(params)
In [2]:
Copied!
figure_path = Path("figures")
figure_path = Path("figures")
In [3]:
Copied!
gamma_catalog = pd.read_csv("gamma_catalog.csv", parse_dates=["time"])
gamma_catalog = pd.read_csv("gamma_catalog.csv", parse_dates=["time"])
In [4]:
Copied!
standard_catalog = pd.read_csv("quakeflow/standard_catalog.csv", parse_dates=["time"])
standard_catalog["depth_km"] = standard_catalog["depth(m)"] / 1e3
standard_catalog_origin = standard_catalog.copy()
standard_catalog = standard_catalog[(standard_catalog["time"] > gamma_catalog["time"].min()) & (standard_catalog["time"] < gamma_catalog["time"].max())]
standard_catalog = pd.read_csv("quakeflow/standard_catalog.csv", parse_dates=["time"])
standard_catalog["depth_km"] = standard_catalog["depth(m)"] / 1e3
standard_catalog_origin = standard_catalog.copy()
standard_catalog = standard_catalog[(standard_catalog["time"] > gamma_catalog["time"].min()) & (standard_catalog["time"] < gamma_catalog["time"].max())]
In [12]:
Copied!
print(len(standard_catalog_origin))
print(len(standard_catalog_origin))
9995
In [ ]:
Copied!
In [5]:
Copied!
seismic_stations = pd.read_json("quakeflow/stations.json", orient="index")
das_loc = pd.read_csv("das_info.csv")
seismic_stations = pd.read_json("quakeflow/stations.json", orient="index")
das_loc = pd.read_csv("das_info.csv")
In [6]:
Copied!
fig, ax = plt.subplots(1, 1, figsize=(4 * 2.5, 3 * 2.5))
ls = LightSource(azdeg=0, altdeg=45)
min_longitude = seismic_stations["longitude"].min() - 0.1
max_longitude = seismic_stations["longitude"].max() + 0.1
min_latitude = seismic_stations["latitude"].min() - 0.1
max_latitude = seismic_stations["latitude"].max() + 0.1
region = [min_longitude, max_longitude, min_latitude, max_latitude]
topo = pygmt.datasets.load_earth_relief(resolution="15s", region=region).to_numpy() / 1e3 # km
x = np.linspace(min_longitude, max_longitude, topo.shape[1])
y = np.linspace(min_latitude, max_latitude, topo.shape[0])
dx, dy = 1, 1
xgrid, ygrid = np.meshgrid(x, y)
ax.pcolormesh(
xgrid,
ygrid,
ls.hillshade(topo, vert_exag=10, dx=dx, dy=dy),
vmin=-1,
shading="gouraud",
cmap="gray",
alpha=1.0,
antialiased=True,
rasterized=True,
)
ax.scatter(das_loc["longitude"], das_loc["latitude"], c="r", s=1, marker=".", label="DAS", rasterized=True)
ax.scatter(seismic_stations["longitude"], seismic_stations["latitude"], c="b", s=30, label="Seismometer", marker="^")
xlim = ax.get_xlim()
ylim = ax.get_ylim()
ax.set_aspect(1.0/np.cos(np.deg2rad(ylim[0])))
plt.savefig(figure_path / "station_map.png")
plt.savefig(figure_path / "station_map.pdf")
fig, ax = plt.subplots(1, 1, figsize=(4 * 2.5, 3 * 2.5))
ls = LightSource(azdeg=0, altdeg=45)
min_longitude = seismic_stations["longitude"].min() - 0.1
max_longitude = seismic_stations["longitude"].max() + 0.1
min_latitude = seismic_stations["latitude"].min() - 0.1
max_latitude = seismic_stations["latitude"].max() + 0.1
region = [min_longitude, max_longitude, min_latitude, max_latitude]
topo = pygmt.datasets.load_earth_relief(resolution="15s", region=region).to_numpy() / 1e3 # km
x = np.linspace(min_longitude, max_longitude, topo.shape[1])
y = np.linspace(min_latitude, max_latitude, topo.shape[0])
dx, dy = 1, 1
xgrid, ygrid = np.meshgrid(x, y)
ax.pcolormesh(
xgrid,
ygrid,
ls.hillshade(topo, vert_exag=10, dx=dx, dy=dy),
vmin=-1,
shading="gouraud",
cmap="gray",
alpha=1.0,
antialiased=True,
rasterized=True,
)
ax.scatter(das_loc["longitude"], das_loc["latitude"], c="r", s=1, marker=".", label="DAS", rasterized=True)
ax.scatter(seismic_stations["longitude"], seismic_stations["latitude"], c="b", s=30, label="Seismometer", marker="^")
xlim = ax.get_xlim()
ylim = ax.get_ylim()
ax.set_aspect(1.0/np.cos(np.deg2rad(ylim[0])))
plt.savefig(figure_path / "station_map.png")
plt.savefig(figure_path / "station_map.pdf")
pygmt-session [WARNING]: Remote dataset given to a data processing module but no registration was specified - default to gridline registration (if available) grdcut [WARNING]: (w - x_min) must equal (NX + eps) * x_inc), where NX is an integer and |eps| <= 0.0001. grdcut [WARNING]: w reset from -119.54162 to -119.541666667 grdcut [WARNING]: (e - x_min) must equal (NX + eps) * x_inc), where NX is an integer and |eps| <= 0.0001. grdcut [WARNING]: e reset from -118.086798 to -118.083333333 grdcut [WARNING]: (s - y_min) must equal (NY + eps) * y_inc), where NY is an integer and |eps| <= 0.0001. grdcut [WARNING]: s reset from 36.84966 to 36.8458333333 grdcut [WARNING]: (n - y_min) must equal (NY + eps) * y_inc), where NY is an integer and |eps| <= 0.0001. grdcut [WARNING]: n reset from 38.4491 to 38.45
In [7]:
Copied!
hypoDD_ct_catalog = pd.read_csv(
"quakeflow/hypoDD_ct_catalog.txt",
sep="\s+",
names=[
"ID",
"LAT",
"LON",
"DEPTH",
"X",
"Y",
"Z",
"EX",
"EY",
"EZ",
"YR",
"MO",
"DY",
"HR",
"MI",
"SC",
"MAG",
"NCCP",
"NCCS",
"NCTP",
"NCTS",
"RCC",
"RCT",
"CID",
],
)
hypoDD_ct_catalog["time"] = hypoDD_ct_catalog.apply(
lambda x: f'{x["YR"]:04.0f}-{x["MO"]:02.0f}-{x["DY"]:02.0f}T{x["HR"]:02.0f}:{x["MI"]:02.0f}:{min(x["SC"], 59.999):05.3f}',
axis=1,
)
hypoDD_ct_catalog["time"] = hypoDD_ct_catalog["time"].apply(lambda x: datetime.strptime(x, "%Y-%m-%dT%H:%M:%S.%f"))
hypoDD_ct_catalog["depth_km"] = hypoDD_ct_catalog["DEPTH"]
hypoDD_ct_catalog["longitude"] = hypoDD_ct_catalog["LON"]
hypoDD_ct_catalog["latitude"] = hypoDD_ct_catalog["LAT"]
hypoDD_ct_catalog = pd.read_csv(
"quakeflow/hypoDD_ct_catalog.txt",
sep="\s+",
names=[
"ID",
"LAT",
"LON",
"DEPTH",
"X",
"Y",
"Z",
"EX",
"EY",
"EZ",
"YR",
"MO",
"DY",
"HR",
"MI",
"SC",
"MAG",
"NCCP",
"NCCS",
"NCTP",
"NCTS",
"RCC",
"RCT",
"CID",
],
)
hypoDD_ct_catalog["time"] = hypoDD_ct_catalog.apply(
lambda x: f'{x["YR"]:04.0f}-{x["MO"]:02.0f}-{x["DY"]:02.0f}T{x["HR"]:02.0f}:{x["MI"]:02.0f}:{min(x["SC"], 59.999):05.3f}',
axis=1,
)
hypoDD_ct_catalog["time"] = hypoDD_ct_catalog["time"].apply(lambda x: datetime.strptime(x, "%Y-%m-%dT%H:%M:%S.%f"))
hypoDD_ct_catalog["depth_km"] = hypoDD_ct_catalog["DEPTH"]
hypoDD_ct_catalog["longitude"] = hypoDD_ct_catalog["LON"]
hypoDD_ct_catalog["latitude"] = hypoDD_ct_catalog["LAT"]
In [8]:
Copied!
quakeflow_catalog_origin = hypoDD_ct_catalog.copy()
quakeflow_catalog_origin = hypoDD_ct_catalog.copy()
In [9]:
Copied!
# hypoDD_ct_catalog2 = hypoDD_ct_catalog[(hypoDD_ct_catalog["longitude"]>min_longitude) & (hypoDD_ct_catalog["longitude"]<max_longitude) & (hypoDD_ct_catalog["latitude"]>min_latitude) & (hypoDD_ct_catalog["latitude"]<max_latitude) & (hypoDD_ct_catalog["depth_km"] > min_depth) & (hypoDD_ct_catalog["depth_km"] < max_depth)]
# plot_location(hypoDD_ct_catalog2, color="k", size=0.3, fname="figures/hypodd_catalog0", xlim=xlim2, ylim=ylim2, zlim=zlim2, stations=das_loc)
# hypoDD_ct_catalog2 = hypoDD_ct_catalog[(hypoDD_ct_catalog["longitude"]>min_longitude) & (hypoDD_ct_catalog["longitude"]min_latitude) & (hypoDD_ct_catalog["latitude"] min_depth) & (hypoDD_ct_catalog["depth_km"] < max_depth)]
# plot_location(hypoDD_ct_catalog2, color="k", size=0.3, fname="figures/hypodd_catalog0", xlim=xlim2, ylim=ylim2, zlim=zlim2, stations=das_loc)
In [10]:
Copied!
hypoDD_ct_catalog = hypoDD_ct_catalog[(hypoDD_ct_catalog["time"] > gamma_catalog["time"].min()) & (hypoDD_ct_catalog["time"] < gamma_catalog["time"].max())]
quakeflow_catalog = hypoDD_ct_catalog
hypoDD_ct_catalog = hypoDD_ct_catalog[(hypoDD_ct_catalog["time"] > gamma_catalog["time"].min()) & (hypoDD_ct_catalog["time"] < gamma_catalog["time"].max())]
quakeflow_catalog = hypoDD_ct_catalog
In [14]:
Copied!
print(len(quakeflow_catalog_origin))
print(len(quakeflow_catalog_origin))
23704
In [11]:
Copied!
def plot_location(catalog, color=None, size=10, cmap="viridis", alpha=0.8, fname="catalog", xlim=None, ylim=None, zlim=None, stations=None):
fig, axes = plt.subplots(2, 2, figsize=(10, 10), gridspec_kw={'width_ratios': [3, 1], 'height_ratios': [3, 1]})
if color is None:
color = "k"
axes[0,0].scatter(catalog["longitude"], catalog["latitude"], s=size, c=color, cmap=cmap, alpha=alpha, rasterized=True)
axes[0,0].set_title(f"Number of events: {len(catalog)}")
axes[0,0].autoscale(tight=True)
if xlim is None:
xlim = axes[0,0].get_xlim()
if ylim is None:
ylim = axes[0,0].get_ylim()
if stations is not None:
axes[0, 0].scatter(stations["longitude"], stations["latitude"], s=10, c="b", marker=".", alpha=0.05, rasterized=True)
axes[0,0].set_xlim(xlim)
axes[0,0].set_ylim(ylim)
axes[0,1].scatter(catalog["depth_km"], catalog["latitude"], s=size, c=color, cmap=cmap, alpha=alpha, rasterized=True)
axes[0,1].autoscale(tight=True)
if ylim is not None:
axes[0,1].set_ylim(ylim)
if zlim is not None:
axes[0,1].set_xlim(zlim)
im = axes[1,0].scatter(catalog["longitude"], catalog["depth_km"], s=size, c=color, cmap=cmap, alpha=alpha, rasterized=True)
axes[1,0].autoscale(tight=True)
if xlim is not None:
axes[1,0].set_xlim(xlim)
if zlim is None:
zlim = axes[1,0].get_ylim()
else:
axes[1,0].set_ylim(zlim)
axes[1,0].invert_yaxis()
# fig.colorbar(im, ax=axes[1,1], location="right", anchor=(-1, -1))
axes[1, 1].axis('off')
fig.tight_layout()
fig.savefig(fname + ".pdf", bbox_inches='tight', dpi=300)
fig.savefig(fname + ".png", bbox_inches='tight', dpi=300)
return xlim, ylim, zlim
def plot_location(catalog, color=None, size=10, cmap="viridis", alpha=0.8, fname="catalog", xlim=None, ylim=None, zlim=None, stations=None):
fig, axes = plt.subplots(2, 2, figsize=(10, 10), gridspec_kw={'width_ratios': [3, 1], 'height_ratios': [3, 1]})
if color is None:
color = "k"
axes[0,0].scatter(catalog["longitude"], catalog["latitude"], s=size, c=color, cmap=cmap, alpha=alpha, rasterized=True)
axes[0,0].set_title(f"Number of events: {len(catalog)}")
axes[0,0].autoscale(tight=True)
if xlim is None:
xlim = axes[0,0].get_xlim()
if ylim is None:
ylim = axes[0,0].get_ylim()
if stations is not None:
axes[0, 0].scatter(stations["longitude"], stations["latitude"], s=10, c="b", marker=".", alpha=0.05, rasterized=True)
axes[0,0].set_xlim(xlim)
axes[0,0].set_ylim(ylim)
axes[0,1].scatter(catalog["depth_km"], catalog["latitude"], s=size, c=color, cmap=cmap, alpha=alpha, rasterized=True)
axes[0,1].autoscale(tight=True)
if ylim is not None:
axes[0,1].set_ylim(ylim)
if zlim is not None:
axes[0,1].set_xlim(zlim)
im = axes[1,0].scatter(catalog["longitude"], catalog["depth_km"], s=size, c=color, cmap=cmap, alpha=alpha, rasterized=True)
axes[1,0].autoscale(tight=True)
if xlim is not None:
axes[1,0].set_xlim(xlim)
if zlim is None:
zlim = axes[1,0].get_ylim()
else:
axes[1,0].set_ylim(zlim)
axes[1,0].invert_yaxis()
# fig.colorbar(im, ax=axes[1,1], location="right", anchor=(-1, -1))
axes[1, 1].axis('off')
fig.tight_layout()
fig.savefig(fname + ".pdf", bbox_inches='tight', dpi=300)
fig.savefig(fname + ".png", bbox_inches='tight', dpi=300)
return xlim, ylim, zlim
In [12]:
Copied!
# sns.set_theme()
# matplotlib.rcParams.update(params)
# plt.style.use('seaborn')
# print(f"{len(gamma_catalog) = }")
# fig, ax = plt.subplots(1, 1, figsize=(15, 10))
# vmin = gamma_catalog["time"].min()
# vmax = gamma_catalog["time"].max()
# vmin = datetime(2020, 11, 1)
# vmax = datetime(2021, 2, 1)
# ax.hist(gamma_catalog["time"], range=(vmin, vmax), bins=(vmax-vmin).days, label=f"PhaseNet-DAS: {len(gamma_catalog)}", alpha=1.0)
# ax.hist(hypoDD_ct_catalog["time"], range=(vmin, vmax), bins=(vmax-vmin).days, label=f"QuakeFlow: {len(hypoDD_ct_catalog)}", alpha=0.8)
# ax.hist(standard_catalog["time"], range=(vmin, vmax), bins=(vmax-vmin).days, label=f"Standard: {len(standard_catalog)}", alpha=0.6)
# ax.set_ylabel("Frequency / day")
# # ax.tick_params(axis="x", labelrotation=15)
# # ax.xaxis.set_major_formatter(mdates.DateFormatter('%b %d %Y'))
# ax.legend()
# ax.autoscale(tight=True)
# fig.autofmt_xdate(rotation=15)
# fig.savefig(figure_path / "earthquake_number.pdf", bbox_inches='tight', dpi=300)
# fig.savefig(figure_path / "earthquake_number.png", bbox_inches='tight', dpi=300)
# sns.set_theme()
# matplotlib.rcParams.update(params)
# plt.style.use('seaborn')
# print(f"{len(gamma_catalog) = }")
# fig, ax = plt.subplots(1, 1, figsize=(15, 10))
# vmin = gamma_catalog["time"].min()
# vmax = gamma_catalog["time"].max()
# vmin = datetime(2020, 11, 1)
# vmax = datetime(2021, 2, 1)
# ax.hist(gamma_catalog["time"], range=(vmin, vmax), bins=(vmax-vmin).days, label=f"PhaseNet-DAS: {len(gamma_catalog)}", alpha=1.0)
# ax.hist(hypoDD_ct_catalog["time"], range=(vmin, vmax), bins=(vmax-vmin).days, label=f"QuakeFlow: {len(hypoDD_ct_catalog)}", alpha=0.8)
# ax.hist(standard_catalog["time"], range=(vmin, vmax), bins=(vmax-vmin).days, label=f"Standard: {len(standard_catalog)}", alpha=0.6)
# ax.set_ylabel("Frequency / day")
# # ax.tick_params(axis="x", labelrotation=15)
# # ax.xaxis.set_major_formatter(mdates.DateFormatter('%b %d %Y'))
# ax.legend()
# ax.autoscale(tight=True)
# fig.autofmt_xdate(rotation=15)
# fig.savefig(figure_path / "earthquake_number.pdf", bbox_inches='tight', dpi=300)
# fig.savefig(figure_path / "earthquake_number.png", bbox_inches='tight', dpi=300)
In [13]:
Copied!
xlim, ylim, zlim = plot_location(gamma_catalog, color="k", size=1, fname="figures/gamma_catalog_das", stations=das_loc)
xlim, ylim, zlim = plot_location(gamma_catalog, color="k", size=1, fname="figures/gamma_catalog_das", stations=das_loc)
In [14]:
Copied!
min_longitude = -119.1
max_longitude = -118.6
min_latitude = 37.45
max_latitude = 37.75
min_depth = -0.1
max_depth = 15
gamma_catalog2 = gamma_catalog[(gamma_catalog["longitude"]>min_longitude) & (gamma_catalog["longitude"]<max_longitude) & (gamma_catalog["latitude"]>min_latitude) & (gamma_catalog["latitude"]<max_latitude) & (gamma_catalog["depth_km"] > min_depth) & (gamma_catalog["depth_km"] < max_depth)]
gamma_catalog2 = gamma_catalog2.sort_values("gamma_score")
color = gamma_catalog2["gamma_score"].clip(np.NaN, gamma_catalog2["gamma_score"].std()*2)
size = gamma_catalog2["gamma_score"].clip(np.NaN, gamma_catalog2["gamma_score"].std()*2)/(gamma_catalog2["gamma_score"].std()*2)*10
xlim2, ylim2, zlim2 = plot_location(gamma_catalog2, color=color, size=size, cmap="binary", fname="figures/gamma_catalog_das2", stations=das_loc)
min_longitude = -119.1
max_longitude = -118.6
min_latitude = 37.45
max_latitude = 37.75
min_depth = -0.1
max_depth = 15
gamma_catalog2 = gamma_catalog[(gamma_catalog["longitude"]>min_longitude) & (gamma_catalog["longitude"]min_latitude) & (gamma_catalog["latitude"] min_depth) & (gamma_catalog["depth_km"] < max_depth)]
gamma_catalog2 = gamma_catalog2.sort_values("gamma_score")
color = gamma_catalog2["gamma_score"].clip(np.NaN, gamma_catalog2["gamma_score"].std()*2)
size = gamma_catalog2["gamma_score"].clip(np.NaN, gamma_catalog2["gamma_score"].std()*2)/(gamma_catalog2["gamma_score"].std()*2)*10
xlim2, ylim2, zlim2 = plot_location(gamma_catalog2, color=color, size=size, cmap="binary", fname="figures/gamma_catalog_das2", stations=das_loc)
In [41]:
Copied!
print(gamma_catalog["time"].min())
print(gamma_catalog["time"].max())
print(gamma_catalog["time"].min())
print(gamma_catalog["time"].max())
2020-10-20 20:02:55.672000 2021-04-12 05:41:15.836000
In [30]:
Copied!
# sns.set_theme("origin")
# sns.reset_orig()
# matplotlib.rcParams.update(params)
min_longitude = -119.2
max_longitude = -118.5
min_latitude = 37.35
max_latitude = 37.8
min_depth = -0.1
max_depth = 15
cmap = "viridis"
gamma_catalog2 = gamma_catalog[(gamma_catalog["longitude"]>min_longitude) & (gamma_catalog["longitude"]<max_longitude) & (gamma_catalog["latitude"]>min_latitude) & (gamma_catalog["latitude"]<max_latitude) & (gamma_catalog["depth_km"] > min_depth) & (gamma_catalog["depth_km"] < max_depth)]
gamma_catalog2 = gamma_catalog2.sort_values("gamma_score")
# events_das2 = events_das[(events_das["longitude"]>min_longitude) & (events_das["longitude"]<max_longitude) & (events_das["latitude"]>min_latitude) & (events_das["latitude"]<max_latitude) & (events_das["depth_km"] > min_depth) & (events_das["depth_km"] < max_depth)]
# # events_das2 = events_das2.sort_values("gamma_score")
# events_mammoth_north2 = events_mammoth_north[(events_mammoth_north["longitude"]>min_longitude) & (events_mammoth_north["longitude"]<max_longitude) & (events_mammoth_north["latitude"]>min_latitude) & (events_mammoth_north["latitude"]<max_latitude) & (events_mammoth_north["depth_km"] > min_depth) & (events_mammoth_north["depth_km"] < max_depth)]
fig, axes = plt.subplots(2, 2, figsize=(12, 12 * (max_latitude-min_latitude)/((max_longitude-min_longitude)*np.cos(np.deg2rad(min_latitude)))), gridspec_kw={'width_ratios': [3, 1], 'height_ratios': [3, 1]})
region = [min_longitude, max_longitude, min_latitude, max_latitude]
topo = pygmt.datasets.load_earth_relief(resolution="03s", region=region).to_numpy() / 1e3 # km
x = np.linspace(min_longitude, max_longitude, topo.shape[1])
y = np.linspace(min_latitude, max_latitude, topo.shape[0])
dx, dy = 1, 1
xgrid, ygrid = np.meshgrid(x, y)
axes[0,0].pcolormesh(
xgrid,
ygrid,
ls.hillshade(topo, vert_exag=10, dx=dx, dy=dy),
vmin=-1,
shading="gouraud",
cmap="gray",
alpha=1.0,
antialiased=True,
rasterized=True,
)
color = gamma_catalog2["gamma_score"].clip(np.NaN, gamma_catalog2["gamma_score"].std()*2)
# color = "red"
size = gamma_catalog2["gamma_score"].clip(np.NaN, gamma_catalog2["gamma_score"].std()*2)/(gamma_catalog2["gamma_score"].std()*2)*15
# cmap="Reds"
cmap = "OrRd"
axes[0,0].scatter(gamma_catalog2["longitude"], gamma_catalog2["latitude"], s=size, c=color, cmap=cmap, linewidth=0, label="DAS", rasterized=True)
axes[0,0].autoscale(tight=True)
xlim = axes[0,0].get_xlim()
ylim = axes[0,0].get_ylim()
# zlim = (0, 21)
# axes[0,0].set_title(f"Number of events: {len(events_das2)}")
axes[0, 0].scatter(das_loc["longitude"], das_loc["latitude"], s=2, c="k", alpha=0.5, marker=".", rasterized=True)
# axes[0, 0].axis("scaled")
axes[0, 0].set_aspect(1.0/np.cos(np.deg2rad(min_latitude)))
axes[0,0].set_xlim(xlim)
axes[0,0].set_ylim(ylim)
# axes[0, 0].legend(markerscale=5)
axes[0,0].set_ylabel("Latitude")
axes[0,1].scatter(gamma_catalog2["depth_km"], gamma_catalog2["latitude"], s=size, c=color, cmap=cmap, linewidth=0, rasterized=True)
axes[0,1].autoscale(tight=True)
axes[0,1].set_ylim(ylim)
axes[0,1].set_xlim([0, max_depth])
axes[0,1].set_xlabel("Depth (km)")
axes[1,0].scatter(gamma_catalog2["longitude"], gamma_catalog2["depth_km"], s=size, c=color, cmap=cmap, linewidth=0, rasterized=True)
axes[1,0].autoscale(tight=True)
axes[1,0].set_xlim(xlim)
axes[1,0].set_ylim([0, max_depth])
axes[1,0].invert_yaxis()
axes[1,0].set_ylabel("Depth (km)")
axes[1,0].set_xlabel("Longitude")
axes[1, 1].axis('off')
# axes[1, 1].set_xlim([0,20])
# axes[1, 1].set_ylim([0,20])
# axes[1, 1].invert_yaxis()
# fig.tight_layout()
fig.savefig(figure_path / "gamma_catalog_das.pdf")
fig.savefig(figure_path / "gamma_catalog_das.png")
# sns.set_theme("origin")
# sns.reset_orig()
# matplotlib.rcParams.update(params)
min_longitude = -119.2
max_longitude = -118.5
min_latitude = 37.35
max_latitude = 37.8
min_depth = -0.1
max_depth = 15
cmap = "viridis"
gamma_catalog2 = gamma_catalog[(gamma_catalog["longitude"]>min_longitude) & (gamma_catalog["longitude"]min_latitude) & (gamma_catalog["latitude"] min_depth) & (gamma_catalog["depth_km"] < max_depth)]
gamma_catalog2 = gamma_catalog2.sort_values("gamma_score")
# events_das2 = events_das[(events_das["longitude"]>min_longitude) & (events_das["longitude"]min_latitude) & (events_das["latitude"] min_depth) & (events_das["depth_km"] < max_depth)]
# # events_das2 = events_das2.sort_values("gamma_score")
# events_mammoth_north2 = events_mammoth_north[(events_mammoth_north["longitude"]>min_longitude) & (events_mammoth_north["longitude"]min_latitude) & (events_mammoth_north["latitude"] min_depth) & (events_mammoth_north["depth_km"] < max_depth)]
fig, axes = plt.subplots(2, 2, figsize=(12, 12 * (max_latitude-min_latitude)/((max_longitude-min_longitude)*np.cos(np.deg2rad(min_latitude)))), gridspec_kw={'width_ratios': [3, 1], 'height_ratios': [3, 1]})
region = [min_longitude, max_longitude, min_latitude, max_latitude]
topo = pygmt.datasets.load_earth_relief(resolution="03s", region=region).to_numpy() / 1e3 # km
x = np.linspace(min_longitude, max_longitude, topo.shape[1])
y = np.linspace(min_latitude, max_latitude, topo.shape[0])
dx, dy = 1, 1
xgrid, ygrid = np.meshgrid(x, y)
axes[0,0].pcolormesh(
xgrid,
ygrid,
ls.hillshade(topo, vert_exag=10, dx=dx, dy=dy),
vmin=-1,
shading="gouraud",
cmap="gray",
alpha=1.0,
antialiased=True,
rasterized=True,
)
color = gamma_catalog2["gamma_score"].clip(np.NaN, gamma_catalog2["gamma_score"].std()*2)
# color = "red"
size = gamma_catalog2["gamma_score"].clip(np.NaN, gamma_catalog2["gamma_score"].std()*2)/(gamma_catalog2["gamma_score"].std()*2)*15
# cmap="Reds"
cmap = "OrRd"
axes[0,0].scatter(gamma_catalog2["longitude"], gamma_catalog2["latitude"], s=size, c=color, cmap=cmap, linewidth=0, label="DAS", rasterized=True)
axes[0,0].autoscale(tight=True)
xlim = axes[0,0].get_xlim()
ylim = axes[0,0].get_ylim()
# zlim = (0, 21)
# axes[0,0].set_title(f"Number of events: {len(events_das2)}")
axes[0, 0].scatter(das_loc["longitude"], das_loc["latitude"], s=2, c="k", alpha=0.5, marker=".", rasterized=True)
# axes[0, 0].axis("scaled")
axes[0, 0].set_aspect(1.0/np.cos(np.deg2rad(min_latitude)))
axes[0,0].set_xlim(xlim)
axes[0,0].set_ylim(ylim)
# axes[0, 0].legend(markerscale=5)
axes[0,0].set_ylabel("Latitude")
axes[0,1].scatter(gamma_catalog2["depth_km"], gamma_catalog2["latitude"], s=size, c=color, cmap=cmap, linewidth=0, rasterized=True)
axes[0,1].autoscale(tight=True)
axes[0,1].set_ylim(ylim)
axes[0,1].set_xlim([0, max_depth])
axes[0,1].set_xlabel("Depth (km)")
axes[1,0].scatter(gamma_catalog2["longitude"], gamma_catalog2["depth_km"], s=size, c=color, cmap=cmap, linewidth=0, rasterized=True)
axes[1,0].autoscale(tight=True)
axes[1,0].set_xlim(xlim)
axes[1,0].set_ylim([0, max_depth])
axes[1,0].invert_yaxis()
axes[1,0].set_ylabel("Depth (km)")
axes[1,0].set_xlabel("Longitude")
axes[1, 1].axis('off')
# axes[1, 1].set_xlim([0,20])
# axes[1, 1].set_ylim([0,20])
# axes[1, 1].invert_yaxis()
# fig.tight_layout()
fig.savefig(figure_path / "gamma_catalog_das.pdf")
fig.savefig(figure_path / "gamma_catalog_das.png")
pygmt-session [WARNING]: Remote dataset given to a data processing module but no registration was specified - default to gridline registration (if available)
In [40]:
Copied!
# sns.reset_orig()
# matplotlib.rcParams.update(params)
min_longitude = -119.2
max_longitude = -118.5
min_latitude = 37.3
max_latitude = 37.8
min_depth = -0.1
max_depth = 15
cmap = "viridis"
gamma_catalog2 = gamma_catalog[(gamma_catalog["longitude"]>min_longitude) & (gamma_catalog["longitude"]<max_longitude) & (gamma_catalog["latitude"]>min_latitude) & (gamma_catalog["latitude"]<max_latitude) & (gamma_catalog["depth_km"] > min_depth) & (gamma_catalog["depth_km"] < max_depth)]
gamma_catalog2 = gamma_catalog2.sort_values("gamma_score")
fig, axes = plt.subplots(2, 2, figsize=(12, 12 * (max_latitude-min_latitude)/((max_longitude-min_longitude)*np.cos(np.deg2rad(min_latitude)))), gridspec_kw={'width_ratios': [3, 1], 'height_ratios': [3, 1]})
region = [min_longitude, max_longitude, min_latitude, max_latitude]
topo = pygmt.datasets.load_earth_relief(resolution="03s", region=region).to_numpy() / 1e3 # km
x = np.linspace(min_longitude, max_longitude, topo.shape[1])
y = np.linspace(min_latitude, max_latitude, topo.shape[0])
dx, dy = 1, 1
xgrid, ygrid = np.meshgrid(x, y)
axes[0,0].pcolormesh(
xgrid,
ygrid,
ls.hillshade(topo, vert_exag=10, dx=dx, dy=dy),
vmin=-1,
shading="gouraud",
cmap="gray",
alpha=1.0,
antialiased=True,
rasterized=True,
)
std = 2.5
color = gamma_catalog2["gamma_score"].clip(np.NaN, gamma_catalog2["gamma_score"].std()*std)
# color = "red"
size = gamma_catalog2["gamma_score"].clip(np.NaN, gamma_catalog2["gamma_score"].std()*std)/(gamma_catalog2["gamma_score"].std()*std)*5
# cmap="Reds"
# cmap = "OrRd"
# cmap = "gray"
cmap = "binary"
axes[0,0].scatter(gamma_catalog2["longitude"], gamma_catalog2["latitude"], s=size, c=color, cmap=cmap, alpha=1.0, label="PhaseNet-DAS", rasterized=True)
axes[0,0].autoscale(tight=True)
xlim = axes[0,0].get_xlim()
ylim = axes[0,0].get_ylim()
# zlim = (0, 21)
# axes[0,0].set_title(f"Number of events: {len(events_das2)}")
axes[0, 0].scatter(das_loc["longitude"], das_loc["latitude"], s=2, c="b", marker=".", alpha=0.05, rasterized=True)
# axes[0, 0].axis("scaled")
axes[0, 0].set_aspect(1.0/np.cos(np.deg2rad(min_latitude)))
axes[0,0].set_xlim(xlim)
axes[0,0].set_ylim(ylim)
# axes[0, 0].legend(markerscale=5)
axes[0, 0].tick_params(axis="x", labelrotation=15)
axes[0,0].set_ylabel("Latitude")
axes[0,1].scatter(gamma_catalog2["depth_km"], gamma_catalog2["latitude"], s=size, c=color, cmap=cmap, alpha=1.0, rasterized=True)
axes[0,1].autoscale(tight=True)
axes[0,1].set_ylim(ylim)
axes[0,1].set_xlim([0, max_depth])
axes[0,1].set_xlabel("Depth (km)")
axes[1,0].scatter(gamma_catalog2["longitude"], gamma_catalog2["depth_km"], s=size, c=color, cmap=cmap, alpha=1.0, rasterized=True)
axes[1,0].autoscale(tight=True)
axes[1,0].set_xlim(xlim)
axes[1,0].set_ylim([0, max_depth])
axes[1,0].invert_yaxis()
axes[1,0].set_ylabel("Depth (km)")
axes[1,0].tick_params(axis="x", labelrotation=15)
axes[1,0].set_xlabel("Longitude")
axes[1, 1].axis('off')
# axes[1, 1].set_xlim([0,20])
# axes[1, 1].set_ylim([0,20])
# axes[1, 1].invert_yaxis()
# fig.tight_layout()
fig.savefig(figure_path / "gamma_catalog_das2.pdf")
fig.savefig(figure_path / "gamma_catalog_das2.png")
# sns.reset_orig()
# matplotlib.rcParams.update(params)
min_longitude = -119.2
max_longitude = -118.5
min_latitude = 37.3
max_latitude = 37.8
min_depth = -0.1
max_depth = 15
cmap = "viridis"
gamma_catalog2 = gamma_catalog[(gamma_catalog["longitude"]>min_longitude) & (gamma_catalog["longitude"]min_latitude) & (gamma_catalog["latitude"] min_depth) & (gamma_catalog["depth_km"] < max_depth)]
gamma_catalog2 = gamma_catalog2.sort_values("gamma_score")
fig, axes = plt.subplots(2, 2, figsize=(12, 12 * (max_latitude-min_latitude)/((max_longitude-min_longitude)*np.cos(np.deg2rad(min_latitude)))), gridspec_kw={'width_ratios': [3, 1], 'height_ratios': [3, 1]})
region = [min_longitude, max_longitude, min_latitude, max_latitude]
topo = pygmt.datasets.load_earth_relief(resolution="03s", region=region).to_numpy() / 1e3 # km
x = np.linspace(min_longitude, max_longitude, topo.shape[1])
y = np.linspace(min_latitude, max_latitude, topo.shape[0])
dx, dy = 1, 1
xgrid, ygrid = np.meshgrid(x, y)
axes[0,0].pcolormesh(
xgrid,
ygrid,
ls.hillshade(topo, vert_exag=10, dx=dx, dy=dy),
vmin=-1,
shading="gouraud",
cmap="gray",
alpha=1.0,
antialiased=True,
rasterized=True,
)
std = 2.5
color = gamma_catalog2["gamma_score"].clip(np.NaN, gamma_catalog2["gamma_score"].std()*std)
# color = "red"
size = gamma_catalog2["gamma_score"].clip(np.NaN, gamma_catalog2["gamma_score"].std()*std)/(gamma_catalog2["gamma_score"].std()*std)*5
# cmap="Reds"
# cmap = "OrRd"
# cmap = "gray"
cmap = "binary"
axes[0,0].scatter(gamma_catalog2["longitude"], gamma_catalog2["latitude"], s=size, c=color, cmap=cmap, alpha=1.0, label="PhaseNet-DAS", rasterized=True)
axes[0,0].autoscale(tight=True)
xlim = axes[0,0].get_xlim()
ylim = axes[0,0].get_ylim()
# zlim = (0, 21)
# axes[0,0].set_title(f"Number of events: {len(events_das2)}")
axes[0, 0].scatter(das_loc["longitude"], das_loc["latitude"], s=2, c="b", marker=".", alpha=0.05, rasterized=True)
# axes[0, 0].axis("scaled")
axes[0, 0].set_aspect(1.0/np.cos(np.deg2rad(min_latitude)))
axes[0,0].set_xlim(xlim)
axes[0,0].set_ylim(ylim)
# axes[0, 0].legend(markerscale=5)
axes[0, 0].tick_params(axis="x", labelrotation=15)
axes[0,0].set_ylabel("Latitude")
axes[0,1].scatter(gamma_catalog2["depth_km"], gamma_catalog2["latitude"], s=size, c=color, cmap=cmap, alpha=1.0, rasterized=True)
axes[0,1].autoscale(tight=True)
axes[0,1].set_ylim(ylim)
axes[0,1].set_xlim([0, max_depth])
axes[0,1].set_xlabel("Depth (km)")
axes[1,0].scatter(gamma_catalog2["longitude"], gamma_catalog2["depth_km"], s=size, c=color, cmap=cmap, alpha=1.0, rasterized=True)
axes[1,0].autoscale(tight=True)
axes[1,0].set_xlim(xlim)
axes[1,0].set_ylim([0, max_depth])
axes[1,0].invert_yaxis()
axes[1,0].set_ylabel("Depth (km)")
axes[1,0].tick_params(axis="x", labelrotation=15)
axes[1,0].set_xlabel("Longitude")
axes[1, 1].axis('off')
# axes[1, 1].set_xlim([0,20])
# axes[1, 1].set_ylim([0,20])
# axes[1, 1].invert_yaxis()
# fig.tight_layout()
fig.savefig(figure_path / "gamma_catalog_das2.pdf")
fig.savefig(figure_path / "gamma_catalog_das2.png")
pygmt-session [WARNING]: Remote dataset given to a data processing module but no registration was specified - default to gridline registration (if available)
In [38]:
Copied!
# sns.reset_orig()
# matplotlib.rcParams.update(params)
min_longitude = -119.2
max_longitude = -118.5
min_latitude = 37.3
max_latitude = 37.8
min_depth = -0.1
max_depth = 15
cmap = "viridis"
quakeflow_catalog2 = quakeflow_catalog[(quakeflow_catalog["longitude"]>min_longitude) & (quakeflow_catalog["longitude"]<max_longitude) & (quakeflow_catalog["latitude"]>min_latitude) & (quakeflow_catalog["latitude"]<max_latitude) & (quakeflow_catalog["depth_km"] > min_depth) & (quakeflow_catalog["depth_km"] < max_depth)]
# quakeflow_catalog2 = quakeflow_catalog2.sort_values("gamma_score")
fig, axes = plt.subplots(2, 2, figsize=(12, 12 * (max_latitude-min_latitude)/((max_longitude-min_longitude)*np.cos(np.deg2rad(min_latitude)))), gridspec_kw={'width_ratios': [3, 1], 'height_ratios': [3, 1]})
region = [min_longitude, max_longitude, min_latitude, max_latitude]
topo = pygmt.datasets.load_earth_relief(resolution="03s", region=region).to_numpy() / 1e3 # km
x = np.linspace(min_longitude, max_longitude, topo.shape[1])
y = np.linspace(min_latitude, max_latitude, topo.shape[0])
dx, dy = 1, 1
xgrid, ygrid = np.meshgrid(x, y)
axes[0,0].pcolormesh(
xgrid,
ygrid,
ls.hillshade(topo, vert_exag=10, dx=dx, dy=dy),
vmin=-1,
shading="gouraud",
cmap="gray",
alpha=1.0,
antialiased=True,
rasterized=True,
)
# std = 2.5
# color = gamma_catalog2["gamma_score"].clip(np.NaN, gamma_catalog2["gamma_score"].std()*std)
# color = "red"
# size = gamma_catalog2["gamma_score"].clip(np.NaN, gamma_catalog2["gamma_score"].std()*std)/(gamma_catalog2["gamma_score"].std()*std)*5
# cmap="Reds"
# cmap = "OrRd"
# cmap = "gray"
size = 2.0
color="k"
cmap="binary"
cmap = "binary"
axes[0,0].scatter(quakeflow_catalog2["longitude"], quakeflow_catalog2["latitude"], s=size, c=color, cmap=cmap, alpha=1.0, label="PhaseNet-DAS", rasterized=True)
axes[0,0].autoscale(tight=True)
xlim = axes[0,0].get_xlim()
ylim = axes[0,0].get_ylim()
# zlim = (0, 21)
# axes[0,0].set_title(f"Number of events: {len(events_das2)}")
axes[0, 0].scatter(seismic_stations["longitude"], seismic_stations["latitude"], s=120, c="b", marker="^", alpha=0.4, rasterized=True)
# axes[0, 0].axis("scaled")
axes[0, 0].set_aspect(1.0/np.cos(np.deg2rad(min_latitude)))
axes[0,0].set_xlim(xlim)
axes[0,0].set_ylim(ylim)
# axes[0, 0].legend(markerscale=5)
axes[0, 0].tick_params(axis="x", labelrotation=15)
axes[0,0].set_ylabel("Latitude")
axes[0,1].scatter(quakeflow_catalog2["depth_km"], quakeflow_catalog2["latitude"], s=size, c=color, cmap=cmap, alpha=1.0, rasterized=True)
axes[0,1].autoscale(tight=True)
axes[0,1].set_ylim(ylim)
axes[0,1].set_xlim([0, max_depth])
axes[0,1].set_xlabel("Depth (km)")
axes[1,0].scatter(quakeflow_catalog2["longitude"], quakeflow_catalog2["depth_km"], s=size, c=color, cmap=cmap, alpha=1.0, rasterized=True)
axes[1,0].autoscale(tight=True)
axes[1,0].set_xlim(xlim)
axes[1,0].set_ylim([0, max_depth])
axes[1,0].invert_yaxis()
axes[1,0].set_ylabel("Depth (km)")
axes[1,0].tick_params(axis="x", labelrotation=15)
axes[1,0].set_xlabel("Longitude")
axes[1, 1].axis('off')
# axes[1, 1].set_xlim([0,20])
# axes[1, 1].set_ylim([0,20])
# axes[1, 1].invert_yaxis()
# fig.tight_layout()
fig.savefig(figure_path / "quakeflow_catalog2.pdf")
fig.savefig(figure_path / "quakeflow_catalog2.png")
# sns.reset_orig()
# matplotlib.rcParams.update(params)
min_longitude = -119.2
max_longitude = -118.5
min_latitude = 37.3
max_latitude = 37.8
min_depth = -0.1
max_depth = 15
cmap = "viridis"
quakeflow_catalog2 = quakeflow_catalog[(quakeflow_catalog["longitude"]>min_longitude) & (quakeflow_catalog["longitude"]min_latitude) & (quakeflow_catalog["latitude"] min_depth) & (quakeflow_catalog["depth_km"] < max_depth)]
# quakeflow_catalog2 = quakeflow_catalog2.sort_values("gamma_score")
fig, axes = plt.subplots(2, 2, figsize=(12, 12 * (max_latitude-min_latitude)/((max_longitude-min_longitude)*np.cos(np.deg2rad(min_latitude)))), gridspec_kw={'width_ratios': [3, 1], 'height_ratios': [3, 1]})
region = [min_longitude, max_longitude, min_latitude, max_latitude]
topo = pygmt.datasets.load_earth_relief(resolution="03s", region=region).to_numpy() / 1e3 # km
x = np.linspace(min_longitude, max_longitude, topo.shape[1])
y = np.linspace(min_latitude, max_latitude, topo.shape[0])
dx, dy = 1, 1
xgrid, ygrid = np.meshgrid(x, y)
axes[0,0].pcolormesh(
xgrid,
ygrid,
ls.hillshade(topo, vert_exag=10, dx=dx, dy=dy),
vmin=-1,
shading="gouraud",
cmap="gray",
alpha=1.0,
antialiased=True,
rasterized=True,
)
# std = 2.5
# color = gamma_catalog2["gamma_score"].clip(np.NaN, gamma_catalog2["gamma_score"].std()*std)
# color = "red"
# size = gamma_catalog2["gamma_score"].clip(np.NaN, gamma_catalog2["gamma_score"].std()*std)/(gamma_catalog2["gamma_score"].std()*std)*5
# cmap="Reds"
# cmap = "OrRd"
# cmap = "gray"
size = 2.0
color="k"
cmap="binary"
cmap = "binary"
axes[0,0].scatter(quakeflow_catalog2["longitude"], quakeflow_catalog2["latitude"], s=size, c=color, cmap=cmap, alpha=1.0, label="PhaseNet-DAS", rasterized=True)
axes[0,0].autoscale(tight=True)
xlim = axes[0,0].get_xlim()
ylim = axes[0,0].get_ylim()
# zlim = (0, 21)
# axes[0,0].set_title(f"Number of events: {len(events_das2)}")
axes[0, 0].scatter(seismic_stations["longitude"], seismic_stations["latitude"], s=120, c="b", marker="^", alpha=0.4, rasterized=True)
# axes[0, 0].axis("scaled")
axes[0, 0].set_aspect(1.0/np.cos(np.deg2rad(min_latitude)))
axes[0,0].set_xlim(xlim)
axes[0,0].set_ylim(ylim)
# axes[0, 0].legend(markerscale=5)
axes[0, 0].tick_params(axis="x", labelrotation=15)
axes[0,0].set_ylabel("Latitude")
axes[0,1].scatter(quakeflow_catalog2["depth_km"], quakeflow_catalog2["latitude"], s=size, c=color, cmap=cmap, alpha=1.0, rasterized=True)
axes[0,1].autoscale(tight=True)
axes[0,1].set_ylim(ylim)
axes[0,1].set_xlim([0, max_depth])
axes[0,1].set_xlabel("Depth (km)")
axes[1,0].scatter(quakeflow_catalog2["longitude"], quakeflow_catalog2["depth_km"], s=size, c=color, cmap=cmap, alpha=1.0, rasterized=True)
axes[1,0].autoscale(tight=True)
axes[1,0].set_xlim(xlim)
axes[1,0].set_ylim([0, max_depth])
axes[1,0].invert_yaxis()
axes[1,0].set_ylabel("Depth (km)")
axes[1,0].tick_params(axis="x", labelrotation=15)
axes[1,0].set_xlabel("Longitude")
axes[1, 1].axis('off')
# axes[1, 1].set_xlim([0,20])
# axes[1, 1].set_ylim([0,20])
# axes[1, 1].invert_yaxis()
# fig.tight_layout()
fig.savefig(figure_path / "quakeflow_catalog2.pdf")
fig.savefig(figure_path / "quakeflow_catalog2.png")
pygmt-session [WARNING]: Remote dataset given to a data processing module but no registration was specified - default to gridline registration (if available)
In [35]:
Copied!
# sns.reset_orig()
# matplotlib.rcParams.update(params)
min_longitude = -119.2
max_longitude = -118.5
min_latitude = 37.3
max_latitude = 37.8
min_depth = -0.1
max_depth = 15
cmap = "viridis"
standard_catalog2 = standard_catalog[(standard_catalog["longitude"]>min_longitude) & (standard_catalog["longitude"]<max_longitude) & (standard_catalog["latitude"]>min_latitude) & (standard_catalog["latitude"]<max_latitude) & (standard_catalog["depth_km"] > min_depth) & (standard_catalog["depth_km"] < max_depth)]
cmap = "viridis_r"
color = standard_catalog2["depth_km"]
fig, axes = plt.subplots(2, 2, figsize=(12, 12 * (max_latitude-min_latitude)/((max_longitude-min_longitude)*np.cos(np.deg2rad(min_latitude)))), gridspec_kw={'width_ratios': [3, 1], 'height_ratios': [3, 1]})
region = [min_longitude, max_longitude, min_latitude, max_latitude]
topo = pygmt.datasets.load_earth_relief(resolution="03s", region=region).to_numpy() / 1e3 # km
x = np.linspace(min_longitude, max_longitude, topo.shape[1])
y = np.linspace(min_latitude, max_latitude, topo.shape[0])
dx, dy = 1, 1
xgrid, ygrid = np.meshgrid(x, y)
axes[0,0].pcolormesh(
xgrid,
ygrid,
ls.hillshade(topo, vert_exag=10, dx=dx, dy=dy),
vmin=-1,
shading="gouraud",
cmap="gray",
alpha=1.0,
antialiased=True,
rasterized=True,
)
# std = 2.5
# color = gamma_catalog2["gamma_score"].clip(np.NaN, gamma_catalog2["gamma_score"].std()*std)
# color = "red"
# size = gamma_catalog2["gamma_score"].clip(np.NaN, gamma_catalog2["gamma_score"].std()*std)/(gamma_catalog2["gamma_score"].std()*std)*5
# cmap="Reds"
# cmap = "OrRd"
# cmap = "gray"
size = 2.0
color="k"
cmap="binary"
cmap = "binary"
axes[0,0].scatter(standard_catalog2["longitude"], standard_catalog2["latitude"], s=size, c=color, cmap=cmap, alpha=1.0, label="PhaseNet-DAS", rasterized=True)
axes[0,0].autoscale(tight=True)
xlim = axes[0,0].get_xlim()
ylim = axes[0,0].get_ylim()
# zlim = (0, 21)
# axes[0,0].set_title(f"Number of events: {len(events_das2)}")
axes[0, 0].scatter(seismic_stations["longitude"], seismic_stations["latitude"], s=120, c="b", marker="^", alpha=0.4, rasterized=True)
# axes[0, 0].axis("scaled")
axes[0, 0].set_aspect(1.0/np.cos(np.deg2rad(min_latitude)))
axes[0,0].set_xlim(xlim)
axes[0,0].set_ylim(ylim)
# axes[0, 0].legend(markerscale=5)
axes[0, 0].tick_params(axis="x", labelrotation=15)
axes[0,0].set_ylabel("Latitude")
axes[0,1].scatter(standard_catalog2["depth_km"], standard_catalog2["latitude"], s=size, c=color, cmap=cmap, alpha=1.0, rasterized=True)
axes[0,1].autoscale(tight=True)
axes[0,1].set_ylim(ylim)
axes[0,1].set_xlim([0, max_depth])
axes[0,1].set_xlabel("Depth (km)")
axes[1,0].scatter(standard_catalog2["longitude"], standard_catalog2["depth_km"], s=size, c=color, cmap=cmap, alpha=1.0, rasterized=True)
axes[1,0].autoscale(tight=True)
axes[1,0].set_xlim(xlim)
axes[1,0].set_ylim([0, max_depth])
axes[1,0].invert_yaxis()
axes[1,0].set_ylabel("Depth (km)")
axes[1,0].tick_params(axis="x", labelrotation=15)
axes[1,0].set_xlabel("Longitude")
axes[1, 1].axis('off')
# axes[1, 1].set_xlim([0,20])
# axes[1, 1].set_ylim([0,20])
# axes[1, 1].invert_yaxis()
# fig.tight_layout()
fig.savefig(figure_path / "standard_catalog2.pdf")
fig.savefig(figure_path / "standard_catalog2.png")
# sns.reset_orig()
# matplotlib.rcParams.update(params)
min_longitude = -119.2
max_longitude = -118.5
min_latitude = 37.3
max_latitude = 37.8
min_depth = -0.1
max_depth = 15
cmap = "viridis"
standard_catalog2 = standard_catalog[(standard_catalog["longitude"]>min_longitude) & (standard_catalog["longitude"]min_latitude) & (standard_catalog["latitude"] min_depth) & (standard_catalog["depth_km"] < max_depth)]
cmap = "viridis_r"
color = standard_catalog2["depth_km"]
fig, axes = plt.subplots(2, 2, figsize=(12, 12 * (max_latitude-min_latitude)/((max_longitude-min_longitude)*np.cos(np.deg2rad(min_latitude)))), gridspec_kw={'width_ratios': [3, 1], 'height_ratios': [3, 1]})
region = [min_longitude, max_longitude, min_latitude, max_latitude]
topo = pygmt.datasets.load_earth_relief(resolution="03s", region=region).to_numpy() / 1e3 # km
x = np.linspace(min_longitude, max_longitude, topo.shape[1])
y = np.linspace(min_latitude, max_latitude, topo.shape[0])
dx, dy = 1, 1
xgrid, ygrid = np.meshgrid(x, y)
axes[0,0].pcolormesh(
xgrid,
ygrid,
ls.hillshade(topo, vert_exag=10, dx=dx, dy=dy),
vmin=-1,
shading="gouraud",
cmap="gray",
alpha=1.0,
antialiased=True,
rasterized=True,
)
# std = 2.5
# color = gamma_catalog2["gamma_score"].clip(np.NaN, gamma_catalog2["gamma_score"].std()*std)
# color = "red"
# size = gamma_catalog2["gamma_score"].clip(np.NaN, gamma_catalog2["gamma_score"].std()*std)/(gamma_catalog2["gamma_score"].std()*std)*5
# cmap="Reds"
# cmap = "OrRd"
# cmap = "gray"
size = 2.0
color="k"
cmap="binary"
cmap = "binary"
axes[0,0].scatter(standard_catalog2["longitude"], standard_catalog2["latitude"], s=size, c=color, cmap=cmap, alpha=1.0, label="PhaseNet-DAS", rasterized=True)
axes[0,0].autoscale(tight=True)
xlim = axes[0,0].get_xlim()
ylim = axes[0,0].get_ylim()
# zlim = (0, 21)
# axes[0,0].set_title(f"Number of events: {len(events_das2)}")
axes[0, 0].scatter(seismic_stations["longitude"], seismic_stations["latitude"], s=120, c="b", marker="^", alpha=0.4, rasterized=True)
# axes[0, 0].axis("scaled")
axes[0, 0].set_aspect(1.0/np.cos(np.deg2rad(min_latitude)))
axes[0,0].set_xlim(xlim)
axes[0,0].set_ylim(ylim)
# axes[0, 0].legend(markerscale=5)
axes[0, 0].tick_params(axis="x", labelrotation=15)
axes[0,0].set_ylabel("Latitude")
axes[0,1].scatter(standard_catalog2["depth_km"], standard_catalog2["latitude"], s=size, c=color, cmap=cmap, alpha=1.0, rasterized=True)
axes[0,1].autoscale(tight=True)
axes[0,1].set_ylim(ylim)
axes[0,1].set_xlim([0, max_depth])
axes[0,1].set_xlabel("Depth (km)")
axes[1,0].scatter(standard_catalog2["longitude"], standard_catalog2["depth_km"], s=size, c=color, cmap=cmap, alpha=1.0, rasterized=True)
axes[1,0].autoscale(tight=True)
axes[1,0].set_xlim(xlim)
axes[1,0].set_ylim([0, max_depth])
axes[1,0].invert_yaxis()
axes[1,0].set_ylabel("Depth (km)")
axes[1,0].tick_params(axis="x", labelrotation=15)
axes[1,0].set_xlabel("Longitude")
axes[1, 1].axis('off')
# axes[1, 1].set_xlim([0,20])
# axes[1, 1].set_ylim([0,20])
# axes[1, 1].invert_yaxis()
# fig.tight_layout()
fig.savefig(figure_path / "standard_catalog2.pdf")
fig.savefig(figure_path / "standard_catalog2.png")
pygmt-session [WARNING]: Remote dataset given to a data processing module but no registration was specified - default to gridline registration (if available)
In [18]:
Copied!
# # sns.reset_orig()
# # matplotlib.rcParams.update(params)
# min_longitude = -119.2
# max_longitude = -118.5
# min_latitude = 37.3
# max_latitude = 37.8
# min_depth = -0.1
# max_depth = 21
# cmap = "viridis"
# standard_catalog2 = standard_catalog[(standard_catalog["longitude"]>min_longitude) & (standard_catalog["longitude"]<max_longitude) & (standard_catalog["latitude"]>min_latitude) & (standard_catalog["latitude"]<max_latitude) & (standard_catalog["depth_km"] > min_depth) & (standard_catalog["depth_km"] < max_depth)]
# cmap = "viridis_r"
# color = standard_catalog2["depth_km"]
# fig, axes = plt.subplots(2, 2, figsize=(12, 12 * (max_latitude-min_latitude)/((max_longitude-min_longitude)*np.cos(np.deg2rad(min_latitude)))), gridspec_kw={'width_ratios': [3, 1], 'height_ratios': [3, 1]})
# region = [min_longitude, max_longitude, min_latitude, max_latitude]
# topo = pygmt.datasets.load_earth_relief(resolution="03s", region=region).to_numpy() / 1e3 # km
# x = np.linspace(min_longitude, max_longitude, topo.shape[1])
# y = np.linspace(min_latitude, max_latitude, topo.shape[0])
# dx, dy = 1, 1
# xgrid, ygrid = np.meshgrid(x, y)
# axes[0,0].pcolormesh(
# xgrid,
# ygrid,
# ls.hillshade(topo, vert_exag=10, dx=dx, dy=dy),
# vmin=-1,
# shading="gouraud",
# cmap="gray",
# alpha=1.0,
# antialiased=True,
# rasterized=True,
# )
# # std = 2.5
# # color = gamma_catalog2["gamma_score"].clip(np.NaN, gamma_catalog2["gamma_score"].std()*std)
# # color = "red"
# # size = gamma_catalog2["gamma_score"].clip(np.NaN, gamma_catalog2["gamma_score"].std()*std)/(gamma_catalog2["gamma_score"].std()*std)*5
# # cmap="Reds"
# # cmap = "OrRd"
# # cmap = "gray"
# size = 2.0
# color="k"
# cmap="binary"
# cmap = "binary"
# axes[0,0].scatter(standard_catalog2["longitude"], standard_catalog2["latitude"], s=size, c=color, cmap=cmap, alpha=1.0, label="PhaseNet-DAS", rasterized=True)
# axes[0,0].autoscale(tight=True)
# xlim = axes[0,0].get_xlim()
# ylim = axes[0,0].get_ylim()
# # zlim = (0, 21)
# # axes[0,0].set_title(f"Number of events: {len(events_das2)}")
# axes[0, 0].scatter(seismic_stations["longitude"], seismic_stations["latitude"], s=120, c="b", marker="^", alpha=0.4, rasterized=True)
# # axes[0, 0].axis("scaled")
# axes[0, 0].set_aspect(1.0/np.cos(np.deg2rad(min_latitude)))
# axes[0,0].set_xlim(xlim)
# axes[0,0].set_ylim(ylim)
# # axes[0, 0].legend(markerscale=5)
# axes[0, 0].tick_params(axis="x", labelrotation=15)
# axes[0,0].set_ylabel("Latitude")
# axes[0,1].scatter(standard_catalog2["depth_km"], standard_catalog2["latitude"], s=size, c=color, cmap=cmap, alpha=1.0, rasterized=True)
# axes[0,1].autoscale(tight=True)
# axes[0,1].set_ylim(ylim)
# axes[0,1].set_xlim([0, 20])
# axes[0,1].set_xlabel("Depth (km)")
# axes[1,0].scatter(standard_catalog2["longitude"], standard_catalog2["depth_km"], s=size, c=color, cmap=cmap, alpha=1.0, rasterized=True)
# axes[1,0].autoscale(tight=True)
# axes[1,0].set_xlim(xlim)
# axes[1,0].set_ylim([0, 20])
# axes[1,0].invert_yaxis()
# axes[1,0].set_ylabel("Depth (km)")
# axes[1,0].tick_params(axis="x", labelrotation=15)
# axes[1,0].set_xlabel("Longitude")
# axes[1, 1].axis('off')
# # axes[1, 1].set_xlim([0,20])
# # axes[1, 1].set_ylim([0,20])
# # axes[1, 1].invert_yaxis()
# # fig.tight_layout()
# fig.savefig(figure_path / "standard_catalog2.pdf")
# fig.savefig(figure_path / "standard_catalog2.png")
# # sns.reset_orig()
# # matplotlib.rcParams.update(params)
# min_longitude = -119.2
# max_longitude = -118.5
# min_latitude = 37.3
# max_latitude = 37.8
# min_depth = -0.1
# max_depth = 21
# cmap = "viridis"
# standard_catalog2 = standard_catalog[(standard_catalog["longitude"]>min_longitude) & (standard_catalog["longitude"]min_latitude) & (standard_catalog["latitude"] min_depth) & (standard_catalog["depth_km"] < max_depth)]
# cmap = "viridis_r"
# color = standard_catalog2["depth_km"]
# fig, axes = plt.subplots(2, 2, figsize=(12, 12 * (max_latitude-min_latitude)/((max_longitude-min_longitude)*np.cos(np.deg2rad(min_latitude)))), gridspec_kw={'width_ratios': [3, 1], 'height_ratios': [3, 1]})
# region = [min_longitude, max_longitude, min_latitude, max_latitude]
# topo = pygmt.datasets.load_earth_relief(resolution="03s", region=region).to_numpy() / 1e3 # km
# x = np.linspace(min_longitude, max_longitude, topo.shape[1])
# y = np.linspace(min_latitude, max_latitude, topo.shape[0])
# dx, dy = 1, 1
# xgrid, ygrid = np.meshgrid(x, y)
# axes[0,0].pcolormesh(
# xgrid,
# ygrid,
# ls.hillshade(topo, vert_exag=10, dx=dx, dy=dy),
# vmin=-1,
# shading="gouraud",
# cmap="gray",
# alpha=1.0,
# antialiased=True,
# rasterized=True,
# )
# # std = 2.5
# # color = gamma_catalog2["gamma_score"].clip(np.NaN, gamma_catalog2["gamma_score"].std()*std)
# # color = "red"
# # size = gamma_catalog2["gamma_score"].clip(np.NaN, gamma_catalog2["gamma_score"].std()*std)/(gamma_catalog2["gamma_score"].std()*std)*5
# # cmap="Reds"
# # cmap = "OrRd"
# # cmap = "gray"
# size = 2.0
# color="k"
# cmap="binary"
# cmap = "binary"
# axes[0,0].scatter(standard_catalog2["longitude"], standard_catalog2["latitude"], s=size, c=color, cmap=cmap, alpha=1.0, label="PhaseNet-DAS", rasterized=True)
# axes[0,0].autoscale(tight=True)
# xlim = axes[0,0].get_xlim()
# ylim = axes[0,0].get_ylim()
# # zlim = (0, 21)
# # axes[0,0].set_title(f"Number of events: {len(events_das2)}")
# axes[0, 0].scatter(seismic_stations["longitude"], seismic_stations["latitude"], s=120, c="b", marker="^", alpha=0.4, rasterized=True)
# # axes[0, 0].axis("scaled")
# axes[0, 0].set_aspect(1.0/np.cos(np.deg2rad(min_latitude)))
# axes[0,0].set_xlim(xlim)
# axes[0,0].set_ylim(ylim)
# # axes[0, 0].legend(markerscale=5)
# axes[0, 0].tick_params(axis="x", labelrotation=15)
# axes[0,0].set_ylabel("Latitude")
# axes[0,1].scatter(standard_catalog2["depth_km"], standard_catalog2["latitude"], s=size, c=color, cmap=cmap, alpha=1.0, rasterized=True)
# axes[0,1].autoscale(tight=True)
# axes[0,1].set_ylim(ylim)
# axes[0,1].set_xlim([0, 20])
# axes[0,1].set_xlabel("Depth (km)")
# axes[1,0].scatter(standard_catalog2["longitude"], standard_catalog2["depth_km"], s=size, c=color, cmap=cmap, alpha=1.0, rasterized=True)
# axes[1,0].autoscale(tight=True)
# axes[1,0].set_xlim(xlim)
# axes[1,0].set_ylim([0, 20])
# axes[1,0].invert_yaxis()
# axes[1,0].set_ylabel("Depth (km)")
# axes[1,0].tick_params(axis="x", labelrotation=15)
# axes[1,0].set_xlabel("Longitude")
# axes[1, 1].axis('off')
# # axes[1, 1].set_xlim([0,20])
# # axes[1, 1].set_ylim([0,20])
# # axes[1, 1].invert_yaxis()
# # fig.tight_layout()
# fig.savefig(figure_path / "standard_catalog2.pdf")
# fig.savefig(figure_path / "standard_catalog2.png")
In [21]:
Copied!
len(quakeflow_catalog_origin2)
len(quakeflow_catalog_origin2)
Out[21]:
17230
In [36]:
Copied!
# sns.reset_orig()
# matplotlib.rcParams.update(params)
min_longitude = -119.2
max_longitude = -118.5
min_latitude = 37.3
max_latitude = 37.8
min_depth = -0.1
max_depth = 15
cmap = "viridis"
quakeflow_catalog_origin2 = quakeflow_catalog_origin[(quakeflow_catalog_origin["longitude"]>min_longitude) & (quakeflow_catalog_origin["longitude"]<max_longitude) & (quakeflow_catalog_origin["latitude"]>min_latitude) & (quakeflow_catalog_origin["latitude"]<max_latitude) & (quakeflow_catalog_origin["depth_km"] > min_depth) & (quakeflow_catalog_origin["depth_km"] < max_depth)]
# quakeflow_catalog2 = quakeflow_catalog2.sort_values("gamma_score")
fig, axes = plt.subplots(2, 2, figsize=(12, 12 * (max_latitude-min_latitude)/((max_longitude-min_longitude)*np.cos(np.deg2rad(min_latitude)))), gridspec_kw={'width_ratios': [3, 1], 'height_ratios': [3, 1]})
region = [min_longitude, max_longitude, min_latitude, max_latitude]
topo = pygmt.datasets.load_earth_relief(resolution="03s", region=region).to_numpy() / 1e3 # km
x = np.linspace(min_longitude, max_longitude, topo.shape[1])
y = np.linspace(min_latitude, max_latitude, topo.shape[0])
dx, dy = 1, 1
xgrid, ygrid = np.meshgrid(x, y)
axes[0,0].pcolormesh(
xgrid,
ygrid,
ls.hillshade(topo, vert_exag=10, dx=dx, dy=dy),
vmin=-1,
shading="gouraud",
cmap="gray",
alpha=1.0,
antialiased=True,
rasterized=True,
)
# std = 2.5
# color = gamma_catalog2["gamma_score"].clip(np.NaN, gamma_catalog2["gamma_score"].std()*std)
# color = "red"
# size = gamma_catalog2["gamma_score"].clip(np.NaN, gamma_catalog2["gamma_score"].std()*std)/(gamma_catalog2["gamma_score"].std()*std)*5
# cmap="Reds"
# cmap = "OrRd"
# cmap = "gray"
size = 0.2
color="k"
cmap="binary"
cmap = "binary"
cmap = "viridis_r"
color = quakeflow_catalog_origin2["depth_km"].clip(np.NaN, 10)
axes[0,0].scatter(quakeflow_catalog_origin2["longitude"], quakeflow_catalog_origin2["latitude"], s=size, c=color, cmap=cmap, alpha=1.0, label="PhaseNet-DAS", rasterized=True)
axes[0,0].autoscale(tight=True)
xlim = axes[0,0].get_xlim()
ylim = axes[0,0].get_ylim()
# zlim = (0, 21)
# axes[0,0].set_title(f"Number of events: {len(events_das2)}")
axes[0, 0].scatter(seismic_stations["longitude"], seismic_stations["latitude"], s=120, c="b", marker="^", alpha=0.4, rasterized=True)
# axes[0, 0].axis("scaled")
axes[0, 0].set_aspect(1.0/np.cos(np.deg2rad(min_latitude)))
axes[0,0].set_xlim(xlim)
axes[0,0].set_ylim(ylim)
# axes[0, 0].legend(markerscale=5)
axes[0, 0].tick_params(axis="x", labelrotation=15)
axes[0,0].set_ylabel("Latitude")
axes[0,1].scatter(quakeflow_catalog_origin2["depth_km"], quakeflow_catalog_origin2["latitude"], s=size, c=color, cmap=cmap, alpha=1.0, rasterized=True)
axes[0,1].autoscale(tight=True)
axes[0,1].set_ylim(ylim)
axes[0,1].set_xlim([0, max_depth])
axes[0,1].set_xlabel("Depth (km)")
axes[1,0].scatter(quakeflow_catalog_origin2["longitude"], quakeflow_catalog_origin2["depth_km"], s=size, c=color, cmap=cmap, alpha=1.0, rasterized=True)
axes[1,0].autoscale(tight=True)
axes[1,0].set_xlim(xlim)
axes[1,0].set_ylim([0, max_depth])
axes[1,0].invert_yaxis()
axes[1,0].set_ylabel("Depth (km)")
axes[1,0].tick_params(axis="x", labelrotation=15)
axes[1,0].set_xlabel("Longitude")
axes[1, 1].axis('off')
# axes[1, 1].set_xlim([0,20])
# axes[1, 1].set_ylim([0,20])
# axes[1, 1].invert_yaxis()
# fig.tight_layout()
fig.savefig(figure_path / "quakeflow_catalog_origin2.pdf")
fig.savefig(figure_path / "quakeflow_catalog_origin2.png")
# sns.reset_orig()
# matplotlib.rcParams.update(params)
min_longitude = -119.2
max_longitude = -118.5
min_latitude = 37.3
max_latitude = 37.8
min_depth = -0.1
max_depth = 15
cmap = "viridis"
quakeflow_catalog_origin2 = quakeflow_catalog_origin[(quakeflow_catalog_origin["longitude"]>min_longitude) & (quakeflow_catalog_origin["longitude"]min_latitude) & (quakeflow_catalog_origin["latitude"] min_depth) & (quakeflow_catalog_origin["depth_km"] < max_depth)]
# quakeflow_catalog2 = quakeflow_catalog2.sort_values("gamma_score")
fig, axes = plt.subplots(2, 2, figsize=(12, 12 * (max_latitude-min_latitude)/((max_longitude-min_longitude)*np.cos(np.deg2rad(min_latitude)))), gridspec_kw={'width_ratios': [3, 1], 'height_ratios': [3, 1]})
region = [min_longitude, max_longitude, min_latitude, max_latitude]
topo = pygmt.datasets.load_earth_relief(resolution="03s", region=region).to_numpy() / 1e3 # km
x = np.linspace(min_longitude, max_longitude, topo.shape[1])
y = np.linspace(min_latitude, max_latitude, topo.shape[0])
dx, dy = 1, 1
xgrid, ygrid = np.meshgrid(x, y)
axes[0,0].pcolormesh(
xgrid,
ygrid,
ls.hillshade(topo, vert_exag=10, dx=dx, dy=dy),
vmin=-1,
shading="gouraud",
cmap="gray",
alpha=1.0,
antialiased=True,
rasterized=True,
)
# std = 2.5
# color = gamma_catalog2["gamma_score"].clip(np.NaN, gamma_catalog2["gamma_score"].std()*std)
# color = "red"
# size = gamma_catalog2["gamma_score"].clip(np.NaN, gamma_catalog2["gamma_score"].std()*std)/(gamma_catalog2["gamma_score"].std()*std)*5
# cmap="Reds"
# cmap = "OrRd"
# cmap = "gray"
size = 0.2
color="k"
cmap="binary"
cmap = "binary"
cmap = "viridis_r"
color = quakeflow_catalog_origin2["depth_km"].clip(np.NaN, 10)
axes[0,0].scatter(quakeflow_catalog_origin2["longitude"], quakeflow_catalog_origin2["latitude"], s=size, c=color, cmap=cmap, alpha=1.0, label="PhaseNet-DAS", rasterized=True)
axes[0,0].autoscale(tight=True)
xlim = axes[0,0].get_xlim()
ylim = axes[0,0].get_ylim()
# zlim = (0, 21)
# axes[0,0].set_title(f"Number of events: {len(events_das2)}")
axes[0, 0].scatter(seismic_stations["longitude"], seismic_stations["latitude"], s=120, c="b", marker="^", alpha=0.4, rasterized=True)
# axes[0, 0].axis("scaled")
axes[0, 0].set_aspect(1.0/np.cos(np.deg2rad(min_latitude)))
axes[0,0].set_xlim(xlim)
axes[0,0].set_ylim(ylim)
# axes[0, 0].legend(markerscale=5)
axes[0, 0].tick_params(axis="x", labelrotation=15)
axes[0,0].set_ylabel("Latitude")
axes[0,1].scatter(quakeflow_catalog_origin2["depth_km"], quakeflow_catalog_origin2["latitude"], s=size, c=color, cmap=cmap, alpha=1.0, rasterized=True)
axes[0,1].autoscale(tight=True)
axes[0,1].set_ylim(ylim)
axes[0,1].set_xlim([0, max_depth])
axes[0,1].set_xlabel("Depth (km)")
axes[1,0].scatter(quakeflow_catalog_origin2["longitude"], quakeflow_catalog_origin2["depth_km"], s=size, c=color, cmap=cmap, alpha=1.0, rasterized=True)
axes[1,0].autoscale(tight=True)
axes[1,0].set_xlim(xlim)
axes[1,0].set_ylim([0, max_depth])
axes[1,0].invert_yaxis()
axes[1,0].set_ylabel("Depth (km)")
axes[1,0].tick_params(axis="x", labelrotation=15)
axes[1,0].set_xlabel("Longitude")
axes[1, 1].axis('off')
# axes[1, 1].set_xlim([0,20])
# axes[1, 1].set_ylim([0,20])
# axes[1, 1].invert_yaxis()
# fig.tight_layout()
fig.savefig(figure_path / "quakeflow_catalog_origin2.pdf")
fig.savefig(figure_path / "quakeflow_catalog_origin2.png")
pygmt-session [WARNING]: Remote dataset given to a data processing module but no registration was specified - default to gridline registration (if available)
In [37]:
Copied!
# sns.reset_orig()
# matplotlib.rcParams.update(params)
min_longitude = -119.2
max_longitude = -118.5
min_latitude = 37.3
max_latitude = 37.8
min_depth = -0.1
max_depth = 15
cmap = "viridis"
standard_catalog_origin2 = standard_catalog_origin[(standard_catalog_origin["longitude"]>min_longitude) & (standard_catalog_origin["longitude"]<max_longitude) & (standard_catalog_origin["latitude"]>min_latitude) & (standard_catalog_origin["latitude"]<max_latitude) & (standard_catalog_origin["depth_km"] > min_depth) & (standard_catalog_origin["depth_km"] < max_depth)]
fig, axes = plt.subplots(2, 2, figsize=(12, 12 * (max_latitude-min_latitude)/((max_longitude-min_longitude)*np.cos(np.deg2rad(min_latitude)))), gridspec_kw={'width_ratios': [3, 1], 'height_ratios': [3, 1]})
region = [min_longitude, max_longitude, min_latitude, max_latitude]
topo = pygmt.datasets.load_earth_relief(resolution="03s", region=region).to_numpy() / 1e3 # km
x = np.linspace(min_longitude, max_longitude, topo.shape[1])
y = np.linspace(min_latitude, max_latitude, topo.shape[0])
dx, dy = 1, 1
xgrid, ygrid = np.meshgrid(x, y)
axes[0,0].pcolormesh(
xgrid,
ygrid,
ls.hillshade(topo, vert_exag=10, dx=dx, dy=dy),
vmin=-1,
shading="gouraud",
cmap="gray",
alpha=1.0,
antialiased=True,
rasterized=True,
)
# std = 2.5
# color = gamma_catalog2["gamma_score"].clip(np.NaN, gamma_catalog2["gamma_score"].std()*std)
# color = "red"
# size = gamma_catalog2["gamma_score"].clip(np.NaN, gamma_catalog2["gamma_score"].std()*std)/(gamma_catalog2["gamma_score"].std()*std)*5
# cmap="Reds"
# cmap = "OrRd"
# cmap = "gray"
size = 0.2
cmap = "viridis_r"
color = standard_catalog_origin2["depth_km"].clip(np.NaN, 10)
axes[0,0].scatter(standard_catalog_origin2["longitude"], standard_catalog_origin2["latitude"], s=size, c=color, cmap=cmap, alpha=1.0, label="PhaseNet-DAS", rasterized=True)
axes[0,0].autoscale(tight=True)
xlim = axes[0,0].get_xlim()
ylim = axes[0,0].get_ylim()
# zlim = (0, 21)
# axes[0,0].set_title(f"Number of events: {len(events_das2)}")
axes[0, 0].scatter(seismic_stations["longitude"], seismic_stations["latitude"], s=120, c="b", marker="^", alpha=0.4, rasterized=True)
# axes[0, 0].axis("scaled")
axes[0, 0].set_aspect(1.0/np.cos(np.deg2rad(min_latitude)))
axes[0,0].set_xlim(xlim)
axes[0,0].set_ylim(ylim)
# axes[0, 0].legend(markerscale=5)
axes[0, 0].tick_params(axis="x", labelrotation=15)
axes[0,0].set_ylabel("Latitude")
axes[0,1].scatter(standard_catalog_origin2["depth_km"], standard_catalog_origin2["latitude"], s=size, c=color, cmap=cmap, alpha=1.0, rasterized=True)
axes[0,1].autoscale(tight=True)
axes[0,1].set_ylim(ylim)
axes[0,1].set_xlim([0, max_depth])
axes[0,1].set_xlabel("Depth (km)")
axes[1,0].scatter(standard_catalog_origin2["longitude"], standard_catalog_origin2["depth_km"], s=size, c=color, cmap=cmap, alpha=1.0, rasterized=True)
axes[1,0].autoscale(tight=True)
axes[1,0].set_xlim(xlim)
axes[1,0].set_ylim([0, max_depth])
axes[1,0].invert_yaxis()
axes[1,0].set_ylabel("Depth (km)")
axes[1,0].tick_params(axis="x", labelrotation=15)
axes[1,0].set_xlabel("Longitude")
axes[1, 1].axis('off')
# axes[1, 1].set_xlim([0,20])
# axes[1, 1].set_ylim([0,20])
# axes[1, 1].invert_yaxis()
# fig.tight_layout()
fig.savefig(figure_path / "standard_catalog_origin2.pdf")
fig.savefig(figure_path / "standard_catalog_origin2.png")
# sns.reset_orig()
# matplotlib.rcParams.update(params)
min_longitude = -119.2
max_longitude = -118.5
min_latitude = 37.3
max_latitude = 37.8
min_depth = -0.1
max_depth = 15
cmap = "viridis"
standard_catalog_origin2 = standard_catalog_origin[(standard_catalog_origin["longitude"]>min_longitude) & (standard_catalog_origin["longitude"]min_latitude) & (standard_catalog_origin["latitude"] min_depth) & (standard_catalog_origin["depth_km"] < max_depth)]
fig, axes = plt.subplots(2, 2, figsize=(12, 12 * (max_latitude-min_latitude)/((max_longitude-min_longitude)*np.cos(np.deg2rad(min_latitude)))), gridspec_kw={'width_ratios': [3, 1], 'height_ratios': [3, 1]})
region = [min_longitude, max_longitude, min_latitude, max_latitude]
topo = pygmt.datasets.load_earth_relief(resolution="03s", region=region).to_numpy() / 1e3 # km
x = np.linspace(min_longitude, max_longitude, topo.shape[1])
y = np.linspace(min_latitude, max_latitude, topo.shape[0])
dx, dy = 1, 1
xgrid, ygrid = np.meshgrid(x, y)
axes[0,0].pcolormesh(
xgrid,
ygrid,
ls.hillshade(topo, vert_exag=10, dx=dx, dy=dy),
vmin=-1,
shading="gouraud",
cmap="gray",
alpha=1.0,
antialiased=True,
rasterized=True,
)
# std = 2.5
# color = gamma_catalog2["gamma_score"].clip(np.NaN, gamma_catalog2["gamma_score"].std()*std)
# color = "red"
# size = gamma_catalog2["gamma_score"].clip(np.NaN, gamma_catalog2["gamma_score"].std()*std)/(gamma_catalog2["gamma_score"].std()*std)*5
# cmap="Reds"
# cmap = "OrRd"
# cmap = "gray"
size = 0.2
cmap = "viridis_r"
color = standard_catalog_origin2["depth_km"].clip(np.NaN, 10)
axes[0,0].scatter(standard_catalog_origin2["longitude"], standard_catalog_origin2["latitude"], s=size, c=color, cmap=cmap, alpha=1.0, label="PhaseNet-DAS", rasterized=True)
axes[0,0].autoscale(tight=True)
xlim = axes[0,0].get_xlim()
ylim = axes[0,0].get_ylim()
# zlim = (0, 21)
# axes[0,0].set_title(f"Number of events: {len(events_das2)}")
axes[0, 0].scatter(seismic_stations["longitude"], seismic_stations["latitude"], s=120, c="b", marker="^", alpha=0.4, rasterized=True)
# axes[0, 0].axis("scaled")
axes[0, 0].set_aspect(1.0/np.cos(np.deg2rad(min_latitude)))
axes[0,0].set_xlim(xlim)
axes[0,0].set_ylim(ylim)
# axes[0, 0].legend(markerscale=5)
axes[0, 0].tick_params(axis="x", labelrotation=15)
axes[0,0].set_ylabel("Latitude")
axes[0,1].scatter(standard_catalog_origin2["depth_km"], standard_catalog_origin2["latitude"], s=size, c=color, cmap=cmap, alpha=1.0, rasterized=True)
axes[0,1].autoscale(tight=True)
axes[0,1].set_ylim(ylim)
axes[0,1].set_xlim([0, max_depth])
axes[0,1].set_xlabel("Depth (km)")
axes[1,0].scatter(standard_catalog_origin2["longitude"], standard_catalog_origin2["depth_km"], s=size, c=color, cmap=cmap, alpha=1.0, rasterized=True)
axes[1,0].autoscale(tight=True)
axes[1,0].set_xlim(xlim)
axes[1,0].set_ylim([0, max_depth])
axes[1,0].invert_yaxis()
axes[1,0].set_ylabel("Depth (km)")
axes[1,0].tick_params(axis="x", labelrotation=15)
axes[1,0].set_xlabel("Longitude")
axes[1, 1].axis('off')
# axes[1, 1].set_xlim([0,20])
# axes[1, 1].set_ylim([0,20])
# axes[1, 1].invert_yaxis()
# fig.tight_layout()
fig.savefig(figure_path / "standard_catalog_origin2.pdf")
fig.savefig(figure_path / "standard_catalog_origin2.png")
pygmt-session [WARNING]: Remote dataset given to a data processing module but no registration was specified - default to gridline registration (if available)
In [ ]:
Copied!
plot_location(standard_catalog, color="k", size=5, fname="figures/standard_catalog", xlim=xlim, ylim=ylim, zlim=zlim, stations=das_loc)
plot_location(standard_catalog, color="k", size=5, fname="figures/standard_catalog", xlim=xlim, ylim=ylim, zlim=zlim, stations=das_loc)
In [ ]:
Copied!
standard_catalog2 = standard_catalog[(standard_catalog["longitude"]>min_longitude) & (standard_catalog["longitude"]<max_longitude) & (standard_catalog["latitude"]>min_latitude) & (standard_catalog["latitude"]<max_latitude) & (standard_catalog["depth_km"] > min_depth) & (standard_catalog["depth_km"] < max_depth)]
plot_location(standard_catalog2, color="k", size=10, fname="figures/standard_catalog2", xlim=xlim2, ylim=ylim2, zlim=zlim2)
standard_catalog2 = standard_catalog[(standard_catalog["longitude"]>min_longitude) & (standard_catalog["longitude"]min_latitude) & (standard_catalog["latitude"] min_depth) & (standard_catalog["depth_km"] < max_depth)]
plot_location(standard_catalog2, color="k", size=10, fname="figures/standard_catalog2", xlim=xlim2, ylim=ylim2, zlim=zlim2)
In [ ]:
Copied!
plot_location(hypoDD_ct_catalog, color="k", size=5, fname="figures/hypodd_catalog", xlim=xlim, ylim=ylim, zlim=zlim, stations=das_loc)
plot_location(hypoDD_ct_catalog, color="k", size=5, fname="figures/hypodd_catalog", xlim=xlim, ylim=ylim, zlim=zlim, stations=das_loc)
In [ ]:
Copied!
hypoDD_ct_catalog2 = hypoDD_ct_catalog[(hypoDD_ct_catalog["longitude"]>min_longitude) & (hypoDD_ct_catalog["longitude"]<max_longitude) & (hypoDD_ct_catalog["latitude"]>min_latitude) & (hypoDD_ct_catalog["latitude"]<max_latitude) & (hypoDD_ct_catalog["depth_km"] > min_depth) & (hypoDD_ct_catalog["depth_km"] < max_depth)]
plot_location(hypoDD_ct_catalog2, color="k", size=10, fname="figures/hypodd_catalog2", xlim=xlim2, ylim=ylim2, zlim=zlim2, stations=das_loc)
hypoDD_ct_catalog2 = hypoDD_ct_catalog[(hypoDD_ct_catalog["longitude"]>min_longitude) & (hypoDD_ct_catalog["longitude"]min_latitude) & (hypoDD_ct_catalog["latitude"] min_depth) & (hypoDD_ct_catalog["depth_km"] < max_depth)]
plot_location(hypoDD_ct_catalog2, color="k", size=10, fname="figures/hypodd_catalog2", xlim=xlim2, ylim=ylim2, zlim=zlim2, stations=das_loc)
In [42]:
Copied!
# picks = pd.read_csv("mammoth_picks.csv", parse_dates=["phase_time"])
gamma_picks = pd.read_csv("gamma_picks.csv", parse_dates=["phase_time"])
# picks = pd.read_csv("mammoth_picks.csv", parse_dates=["phase_time"])
gamma_picks = pd.read_csv("gamma_picks.csv", parse_dates=["phase_time"])
In [72]:
Copied!
gamma_picks = gamma_picks.sort_values("phase_time")
gamma_picks = gamma_picks.sort_values("phase_time")
In [ ]:
Copied!
# fig, axes = plt.subplots(3, 1, figsize=(15, 10), sharex=True, sharey=True)
# axes[0].hist(picks["channel_index"], bins=100)
# axes[0].autoscale(tight=True)
# axes[1].hist(picks[picks["phase_type"] == "P"]["channel_index"], bins=100)
# axes[1].autoscale(tight=True)
# axes[2].hist(picks[picks["phase_type"] == "S"]["channel_index"], bins=100)
# axes[2].autoscale(tight=True)
# fig, axes = plt.subplots(3, 1, figsize=(15, 10), sharex=True, sharey=True)
# axes[0].hist(picks["channel_index"], bins=100)
# axes[0].autoscale(tight=True)
# axes[1].hist(picks[picks["phase_type"] == "P"]["channel_index"], bins=100)
# axes[1].autoscale(tight=True)
# axes[2].hist(picks[picks["phase_type"] == "S"]["channel_index"], bins=100)
# axes[2].autoscale(tight=True)
In [93]:
Copied!
# das_loc["picks"] = das_loc["index"].apply(lambda x: len(picks[picks["channel_index"] == x]))
gamma_picks2 = gamma_picks[gamma_picks["event_id"] != -1]
das_loc["picks"] = das_loc["index"].apply(lambda x: len(gamma_picks2[gamma_picks2["channel_index"] == x]))
# das_loc["picks"] = das_loc["index"].apply(lambda x: len(picks[picks["channel_index"] == x]))
gamma_picks2 = gamma_picks[gamma_picks["event_id"] != -1]
das_loc["picks"] = das_loc["index"].apply(lambda x: len(gamma_picks2[gamma_picks2["channel_index"] == x]))
In [94]:
Copied!
min_latitude, max_latitude = das_loc["latitude"].min()-0.1, das_loc["latitude"].max()+0.1
min_longitude, max_longitude = das_loc["longitude"].min()-0.1, das_loc["longitude"].max()+0.1
# fig, ax = plt.subplots(1, 1, figsize=(10, 10*(max_latitude - min_latitude)/((max_longitude - min_longitude)*np.cos(np.deg2rad(min_latitude)))))
fig, ax = plt.subplots(1, 1, figsize=(10, 10*(max_latitude - min_latitude)/(max_longitude - min_longitude)))
region = [min_longitude, max_longitude, min_latitude, max_latitude]
topo = pygmt.datasets.load_earth_relief(resolution="03s", region=region).to_numpy() / 1e3 # km
x = np.linspace(min_longitude, max_longitude, topo.shape[1])
y = np.linspace(min_latitude, max_latitude, topo.shape[0])
dx, dy = 1, 1
xgrid, ygrid = np.meshgrid(x, y)
ax.pcolormesh(
xgrid,
ygrid,
ls.hillshade(topo, vert_exag=10, dx=dx, dy=dy),
vmin=-1,
shading="gouraud",
cmap="gray",
alpha=1.0,
antialiased=True,
rasterized=True,
)
im = ax.scatter(das_loc["longitude"], das_loc["latitude"], c=das_loc["picks"], label="DAS", vmax=3*np.std(das_loc["picks"]), cmap="viridis")
ax.set_xlabel("Longitude")
ax.set_ylabel("Latitude")
ax.tick_params(axis="x", labelrotation=15)
ax.set_aspect(1.0/np.cos(np.deg2rad(min_latitude)))
clb = fig.colorbar(im, ax=ax)
clb.ax.set_ylabel('Frequency')
# fig.tight_layout()
fig.savefig(figure_path / "pick_frequency_along_das.pdf")
fig.savefig(figure_path / "pick_frequency_along_das.png")
min_latitude, max_latitude = das_loc["latitude"].min()-0.1, das_loc["latitude"].max()+0.1
min_longitude, max_longitude = das_loc["longitude"].min()-0.1, das_loc["longitude"].max()+0.1
# fig, ax = plt.subplots(1, 1, figsize=(10, 10*(max_latitude - min_latitude)/((max_longitude - min_longitude)*np.cos(np.deg2rad(min_latitude)))))
fig, ax = plt.subplots(1, 1, figsize=(10, 10*(max_latitude - min_latitude)/(max_longitude - min_longitude)))
region = [min_longitude, max_longitude, min_latitude, max_latitude]
topo = pygmt.datasets.load_earth_relief(resolution="03s", region=region).to_numpy() / 1e3 # km
x = np.linspace(min_longitude, max_longitude, topo.shape[1])
y = np.linspace(min_latitude, max_latitude, topo.shape[0])
dx, dy = 1, 1
xgrid, ygrid = np.meshgrid(x, y)
ax.pcolormesh(
xgrid,
ygrid,
ls.hillshade(topo, vert_exag=10, dx=dx, dy=dy),
vmin=-1,
shading="gouraud",
cmap="gray",
alpha=1.0,
antialiased=True,
rasterized=True,
)
im = ax.scatter(das_loc["longitude"], das_loc["latitude"], c=das_loc["picks"], label="DAS", vmax=3*np.std(das_loc["picks"]), cmap="viridis")
ax.set_xlabel("Longitude")
ax.set_ylabel("Latitude")
ax.tick_params(axis="x", labelrotation=15)
ax.set_aspect(1.0/np.cos(np.deg2rad(min_latitude)))
clb = fig.colorbar(im, ax=ax)
clb.ax.set_ylabel('Frequency')
# fig.tight_layout()
fig.savefig(figure_path / "pick_frequency_along_das.pdf")
fig.savefig(figure_path / "pick_frequency_along_das.png")
pygmt-session [WARNING]: Remote dataset given to a data processing module but no registration was specified - default to gridline registration (if available) grdcut [WARNING]: (w - x_min) must equal (NX + eps) * x_inc), where NX is an integer and |eps| <= 0.0001. grdcut [WARNING]: w reset from -119.22032423 to -119.220833333 grdcut [WARNING]: (e - x_min) must equal (NX + eps) * x_inc), where NX is an integer and |eps| <= 0.0001. grdcut [WARNING]: e reset from -118.523626 to -118.523333333 grdcut [WARNING]: (s - y_min) must equal (NY + eps) * y_inc), where NY is an integer and |eps| <= 0.0001. grdcut [WARNING]: s reset from 37.4561743054 to 37.4558333333 grdcut [WARNING]: (n - y_min) must equal (NY + eps) * y_inc), where NY is an integer and |eps| <= 0.0001. grdcut [WARNING]: n reset from 38.0619240402 to 38.0625
In [90]:
Copied!
gamma_picks2 = gamma_picks.iloc[:1000_000]
gamma_picks2 = gamma_picks2.sort_values("phase_score")
fig, axes = plt.subplots(2, 1, figsize=(15, 10), sharex=True, sharey=True)
axes[0].scatter(gamma_picks2[gamma_picks2["phase_type"] == "P"]["phase_time"], gamma_picks2[gamma_picks2["phase_type"] == "P"]["channel_index"], c=gamma_picks2[gamma_picks2["phase_type"] == "P"]["phase_score"], s=1, cmap="hot_r", vmin=0.6, vmax=1.0, label="P-picks", rasterized=True)
axes[0].autoscale(tight=True)
axes[0].set_ylim([0, 10_000])
axes[0].set_ylabel("Channel index")
axes[0].legend(loc="lower left")
xlim = axes[0].get_xlim()
axes[1].scatter(gamma_picks2[gamma_picks2["phase_type"] == "S"]["phase_time"], gamma_picks2[gamma_picks2["phase_type"] == "S"]["channel_index"], c=gamma_picks2[gamma_picks2["phase_type"] == "S"]["phase_score"], s=1, cmap="hot_r", vmin=0.6, vmax=1.0, label="S-picks", rasterized=True)
axes[1].autoscale(tight=True)
axes[1].set_ylim([0, 10_000])
axes[1].set_ylabel("Channel index")
axes[1].legend(loc="lower left")
axes[1].set_xlim(xlim)
fig.autofmt_xdate(rotation=15)
fig.tight_layout()
fig.savefig(figure_path / "pick_time_index.png")
fig.savefig(figure_path / "pick_time_index.pdf")
gamma_picks2 = gamma_picks.iloc[:1000_000]
gamma_picks2 = gamma_picks2.sort_values("phase_score")
fig, axes = plt.subplots(2, 1, figsize=(15, 10), sharex=True, sharey=True)
axes[0].scatter(gamma_picks2[gamma_picks2["phase_type"] == "P"]["phase_time"], gamma_picks2[gamma_picks2["phase_type"] == "P"]["channel_index"], c=gamma_picks2[gamma_picks2["phase_type"] == "P"]["phase_score"], s=1, cmap="hot_r", vmin=0.6, vmax=1.0, label="P-picks", rasterized=True)
axes[0].autoscale(tight=True)
axes[0].set_ylim([0, 10_000])
axes[0].set_ylabel("Channel index")
axes[0].legend(loc="lower left")
xlim = axes[0].get_xlim()
axes[1].scatter(gamma_picks2[gamma_picks2["phase_type"] == "S"]["phase_time"], gamma_picks2[gamma_picks2["phase_type"] == "S"]["channel_index"], c=gamma_picks2[gamma_picks2["phase_type"] == "S"]["phase_score"], s=1, cmap="hot_r", vmin=0.6, vmax=1.0, label="S-picks", rasterized=True)
axes[1].autoscale(tight=True)
axes[1].set_ylim([0, 10_000])
axes[1].set_ylabel("Channel index")
axes[1].legend(loc="lower left")
axes[1].set_xlim(xlim)
fig.autofmt_xdate(rotation=15)
fig.tight_layout()
fig.savefig(figure_path / "pick_time_index.png")
fig.savefig(figure_path / "pick_time_index.pdf")
In [98]:
Copied!
gamma_picks2 = gamma_picks.iloc[::100]
gamma_picks2 = gamma_picks2.sort_values("phase_score")
fig, axes = plt.subplots(2, 1, figsize=(15, 10), sharex=True, sharey=True)
axes[0].scatter(gamma_picks2[gamma_picks2["phase_type"] == "P"]["phase_time"], gamma_picks2[gamma_picks2["phase_type"] == "P"]["channel_index"], c=gamma_picks2[gamma_picks2["phase_type"] == "P"]["phase_score"], s=1, cmap="hot_r", vmin=0.6, vmax=1.0, label="P-picks", rasterized=True)
axes[0].autoscale(tight=True)
axes[0].set_ylim([0, 10_000])
axes[0].set_ylabel("Channel index")
axes[0].legend(loc="lower left")
xlim = axes[0].get_xlim()
axes[1].scatter(gamma_picks2[gamma_picks2["phase_type"] == "S"]["phase_time"], gamma_picks2[gamma_picks2["phase_type"] == "S"]["channel_index"], c=gamma_picks2[gamma_picks2["phase_type"] == "S"]["phase_score"], s=1, cmap="hot_r", vmin=0.6, vmax=1.0, label="S-picks", rasterized=True)
axes[1].autoscale(tight=True)
axes[1].set_ylim([0, 10_000])
axes[1].set_ylabel("Channel index")
axes[1].legend(loc="lower left")
axes[1].set_xlim(xlim)
fig.autofmt_xdate(rotation=15)
fig.tight_layout()
fig.savefig(figure_path / "pick_time_index_origin.png")
fig.savefig(figure_path / "pick_time_index_origin.pdf")
gamma_picks2 = gamma_picks.iloc[::100]
gamma_picks2 = gamma_picks2.sort_values("phase_score")
fig, axes = plt.subplots(2, 1, figsize=(15, 10), sharex=True, sharey=True)
axes[0].scatter(gamma_picks2[gamma_picks2["phase_type"] == "P"]["phase_time"], gamma_picks2[gamma_picks2["phase_type"] == "P"]["channel_index"], c=gamma_picks2[gamma_picks2["phase_type"] == "P"]["phase_score"], s=1, cmap="hot_r", vmin=0.6, vmax=1.0, label="P-picks", rasterized=True)
axes[0].autoscale(tight=True)
axes[0].set_ylim([0, 10_000])
axes[0].set_ylabel("Channel index")
axes[0].legend(loc="lower left")
xlim = axes[0].get_xlim()
axes[1].scatter(gamma_picks2[gamma_picks2["phase_type"] == "S"]["phase_time"], gamma_picks2[gamma_picks2["phase_type"] == "S"]["channel_index"], c=gamma_picks2[gamma_picks2["phase_type"] == "S"]["phase_score"], s=1, cmap="hot_r", vmin=0.6, vmax=1.0, label="S-picks", rasterized=True)
axes[1].autoscale(tight=True)
axes[1].set_ylim([0, 10_000])
axes[1].set_ylabel("Channel index")
axes[1].legend(loc="lower left")
axes[1].set_xlim(xlim)
fig.autofmt_xdate(rotation=15)
fig.tight_layout()
fig.savefig(figure_path / "pick_time_index_origin.png")
fig.savefig(figure_path / "pick_time_index_origin.pdf")
In [91]:
Copied!
gamma_picks2 = gamma_picks.iloc[:1000_000]
gamma_picks2 = gamma_picks2.sort_values("phase_score")
gamma_picks2 = gamma_picks2[gamma_picks2["event_id"] != -1]
fig, axes = plt.subplots(2, 1, figsize=(15, 10), sharex=True, sharey=True)
axes[0].scatter(gamma_picks2[gamma_picks2["phase_type"] == "P"]["phase_time"], gamma_picks2[gamma_picks2["phase_type"] == "P"]["channel_index"], c=gamma_picks2[gamma_picks2["phase_type"] == "P"]["phase_score"], s=1, cmap="hot_r", vmin=0.6, vmax=1.0, label="P-picks", rasterized=True)
axes[0].autoscale(tight=True)
axes[0].set_ylim([0, 10_000])
axes[0].set_ylabel("Channel index")
axes[0].legend(loc="lower left")
axes[0].set_xlim(xlim)
axes[1].scatter(gamma_picks2[gamma_picks2["phase_type"] == "S"]["phase_time"], gamma_picks2[gamma_picks2["phase_type"] == "S"]["channel_index"], c=gamma_picks2[gamma_picks2["phase_type"] == "S"]["phase_score"], s=1, cmap="hot_r", vmin=0.6, vmax=1.0, label="S-picks", rasterized=True)
axes[1].autoscale(tight=True)
axes[1].set_ylim([0, 10_000])
axes[1].set_ylabel("Channel index")
axes[1].legend(loc="lower left")
axes[1].set_xlim(xlim)
fig.autofmt_xdate(rotation=15)
fig.tight_layout()
fig.savefig(figure_path / "pick_time_index_gamma.png")
fig.savefig(figure_path / "pick_time_index_gamma.pdf")
gamma_picks2 = gamma_picks.iloc[:1000_000]
gamma_picks2 = gamma_picks2.sort_values("phase_score")
gamma_picks2 = gamma_picks2[gamma_picks2["event_id"] != -1]
fig, axes = plt.subplots(2, 1, figsize=(15, 10), sharex=True, sharey=True)
axes[0].scatter(gamma_picks2[gamma_picks2["phase_type"] == "P"]["phase_time"], gamma_picks2[gamma_picks2["phase_type"] == "P"]["channel_index"], c=gamma_picks2[gamma_picks2["phase_type"] == "P"]["phase_score"], s=1, cmap="hot_r", vmin=0.6, vmax=1.0, label="P-picks", rasterized=True)
axes[0].autoscale(tight=True)
axes[0].set_ylim([0, 10_000])
axes[0].set_ylabel("Channel index")
axes[0].legend(loc="lower left")
axes[0].set_xlim(xlim)
axes[1].scatter(gamma_picks2[gamma_picks2["phase_type"] == "S"]["phase_time"], gamma_picks2[gamma_picks2["phase_type"] == "S"]["channel_index"], c=gamma_picks2[gamma_picks2["phase_type"] == "S"]["phase_score"], s=1, cmap="hot_r", vmin=0.6, vmax=1.0, label="S-picks", rasterized=True)
axes[1].autoscale(tight=True)
axes[1].set_ylim([0, 10_000])
axes[1].set_ylabel("Channel index")
axes[1].legend(loc="lower left")
axes[1].set_xlim(xlim)
fig.autofmt_xdate(rotation=15)
fig.tight_layout()
fig.savefig(figure_path / "pick_time_index_gamma.png")
fig.savefig(figure_path / "pick_time_index_gamma.pdf")
In [129]:
Copied!
gamma_picks2 = gamma_picks.iloc[:1000_000]
gamma_picks2 = gamma_picks2.sort_values("phase_score")
gamma_picks2 = gamma_picks2[gamma_picks2["event_id"] != -1]
fig, axes = plt.subplots(2, 1, figsize=(10, 10), sharex=True, sharey=True)
# tmin = datetime(2020, 10, 22, 16, 30, )
# tmax = datetime(2020, 10, 22, 16, 30, 20)
# tmin = datetime(2020, 10, 28, 1, 30, )
# tmax = datetime(2020, 10, 28, 2, 1, )
tmin = datetime(2020, 10, 28, 1, 42, 40)
tmax = datetime(2020, 10, 28, 1, 42, 55)
axes[0].scatter(gamma_picks2[gamma_picks2["phase_type"] == "P"]["phase_time"], gamma_picks2[gamma_picks2["phase_type"] == "P"]["channel_index"], c=gamma_picks2[gamma_picks2["phase_type"] == "P"]["phase_score"], s=50, cmap="hot_r", vmin=0.6, vmax=1.0, label="P-picks", rasterized=True)
axes[0].autoscale(tight=True)
# axes[0].set_ylim([0, 10_000])
axes[0].set_ylabel("Channel index")
axes[0].legend(loc="lower left")
axes[0].set_xlim([tmin, tmax])
axes[1].scatter(gamma_picks2[gamma_picks2["phase_type"] == "S"]["phase_time"], gamma_picks2[gamma_picks2["phase_type"] == "S"]["channel_index"], c=gamma_picks2[gamma_picks2["phase_type"] == "S"]["phase_score"], s=50, cmap="hot_r", vmin=0.6, vmax=1.0, label="S-picks", rasterized=True)
axes[1].autoscale(tight=True)
# axes[1].set_ylim([0, 10_000])
axes[1].set_ylabel("Channel index")
axes[1].legend(loc="lower left")
axes[1].set_xlim([tmin, tmax])
fig.autofmt_xdate(rotation=15)
fig.tight_layout()
fig.savefig(figure_path / "pick_time_index2.png")
fig.savefig(figure_path / "pick_time_index2.pdf")
gamma_picks2 = gamma_picks.iloc[:1000_000]
gamma_picks2 = gamma_picks2.sort_values("phase_score")
gamma_picks2 = gamma_picks2[gamma_picks2["event_id"] != -1]
fig, axes = plt.subplots(2, 1, figsize=(10, 10), sharex=True, sharey=True)
# tmin = datetime(2020, 10, 22, 16, 30, )
# tmax = datetime(2020, 10, 22, 16, 30, 20)
# tmin = datetime(2020, 10, 28, 1, 30, )
# tmax = datetime(2020, 10, 28, 2, 1, )
tmin = datetime(2020, 10, 28, 1, 42, 40)
tmax = datetime(2020, 10, 28, 1, 42, 55)
axes[0].scatter(gamma_picks2[gamma_picks2["phase_type"] == "P"]["phase_time"], gamma_picks2[gamma_picks2["phase_type"] == "P"]["channel_index"], c=gamma_picks2[gamma_picks2["phase_type"] == "P"]["phase_score"], s=50, cmap="hot_r", vmin=0.6, vmax=1.0, label="P-picks", rasterized=True)
axes[0].autoscale(tight=True)
# axes[0].set_ylim([0, 10_000])
axes[0].set_ylabel("Channel index")
axes[0].legend(loc="lower left")
axes[0].set_xlim([tmin, tmax])
axes[1].scatter(gamma_picks2[gamma_picks2["phase_type"] == "S"]["phase_time"], gamma_picks2[gamma_picks2["phase_type"] == "S"]["channel_index"], c=gamma_picks2[gamma_picks2["phase_type"] == "S"]["phase_score"], s=50, cmap="hot_r", vmin=0.6, vmax=1.0, label="S-picks", rasterized=True)
axes[1].autoscale(tight=True)
# axes[1].set_ylim([0, 10_000])
axes[1].set_ylabel("Channel index")
axes[1].legend(loc="lower left")
axes[1].set_xlim([tmin, tmax])
fig.autofmt_xdate(rotation=15)
fig.tight_layout()
fig.savefig(figure_path / "pick_time_index2.png")
fig.savefig(figure_path / "pick_time_index2.pdf")
In [ ]:
Copied!